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Introduction 

In  this  lecture  notes  we  review  some  recent  results  concerning  the  numerical  solution  of 
nonlinear  collisional  kinetic  equation.  The  most  well-known  example  is  represented  by  the 
Boltzmann  equation  of  rarefied  gas  dynamics  (Cercignani,  1988;  Cercignani  et  ah,  1994). 
Besides  other  classical  examples,  like  the  Landau  equation  of  plasma  physics  (Landau, 
1981),  kinetic  equations  play  an  important  role  in  modelling  granular  gases  (Bobylev  et  al., 
2000),  charged  particles  in  semiconductors  (Markowich  et  al.,  1989),  neutron  transport 
(Jin  et  ah,  2000)  and  quantum  gases  (Escobedo  et  ah,  2003b).  More  recently  applica¬ 
tions  of  kinetic  equations  have  been  considered  for  car  traffic  flows  (Klar  and  Wegener, 
1997),  chemotactical  movements  (Chalub  et  ah,  2004),  tumor  immune  cells  competition 
(Bcllomo  and  Bcllouquid,  2004),  coagulation- fragmentation  processes  (Escobedo  et  ah, 
2003a),  population  dynamics  (Desvillettes  et  ah,  2004),  market  economies  (Cordier  et  ah, 
2005),  supply  chains  (Armbruster  et  ah,  2007),  flocking  dynamics  (Ha  and  Tadmor,  2008) 
and  many  other.  For  a  recent  introduction  to  the  Boltzmann  equation  and  related  kinetic 
equations  we  refer  the  reader  to  Degond  et  ah  (2004);  Villani  (2002),  recent  applications 
to  biology  and  socio-economy  can  be  found  in  Naldi  et  ah  (2010). 

Although  the  scope  of  our  insights  is  wider,  here  we  will  focus  mainly  on  the  classical 
Boltzmann  equation  of  rarefied  gas  dynamics.  This  is  motivated  not  only  by  its  relevance 
for  applications  but  also  because  it  contains  all  major  difficulties  present  in  other  kinetic 
models  and  represents  the  most  challenging  case  for  the  development  of  numerical  schemes. 

Approximate  methods  of  solution  for  the  Boltzmann  equation  have  a  long  history 
tracing  back  to  Hilbert,  Chapmann  and  Enskog  (Cercignani,  1988)  at  the  beginning  of 
the  last  century.  The  mathematical  difficulties  related  to  the  Boltzmann  equation  make 
it  extremely  difficult,  if  not  impossible,  the  determination  of  analytic  solutions  in  most 
physically  relevant  situations.  Only  in  recent  years,  starting  in  the  70s  with  the  pioneering 
works  by  Chorin  (1972)  and  Sod  (1977),  the  problem  has  been  tackled  numerically  with 
particular  care  to  accuracy  and  computational  cost. 

Most  of  the  difficulties  are  due  to  the  multidimensional  structure  of  the  collisional 
integral,  since  the  integration  runs  on  a  highly-dimensional  unflat  manifold.  In  addition 
the  numerical  integration  requires  great  care  since  the  collision  integral  is  at  the  basis 
of  the  macroscopic  properties  of  the  equation.  Further  difficulties  are  represented  by  the 
presence  of  stiffness,  like  the  case  of  small  mean  free  path  (Gabetta  et  al.,  1997)  or  the 
case  of  large  velocities  (Filbet  and  Pareschi,  2003). 

For  such  reasons  realistic  numerical  simulations  are  based  on  Monte-Carlo  techniques. 
The  most  famous  examples  are  the  Direct  Simulation  Monte-Carlo  (DSMC)  methods  by 
Bird  (Bird,  1994)  and  by  Nanbu  (Nanbu,  1980).  These  methods  guarantee  efficiency  and 
preservation  of  the  main  physical  properties.  However,  avoiding  statistical  fluctuations  in 
the  results  becomes  extremely  expensive  in  presence  of  non-stationary  flows  or  close  to 
continuum  regimes. 

Among  deterministic  approximations  one  of  the  most  popular  method  is  represented 
by  the  so  called  Discrete  Velocity  Models  (DVM)  of  the  Boltzmann  equation.  These  meth¬ 
ods  (Martin  et  ah,  1992;  Rogier  and  Schneider,  1994;  Bobylev  et  ah,  1995;  Buet,  1996; 
Panferov  and  Heintz,  2002)  are  based  on  a  Cartesian  grid  in  velocity  and  on  a  discrete 
collision  mechanism  on  the  points  of  the  grid  that  preserves  the  main  physical  proper¬ 
ties.  LInfortunately  DVM  are  not  competitive  with  Monte  Carlo  methods  in  terms  of 
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computational  cost  (typically  0(n^2dv+1^dv),  where  n  is  the  total  number  of  discretization 
parameters  in  velocity  and  dv  is  the  dimension  of  the  velocity  space)  and  their  accu¬ 
racy  seems  to  be  at  most  first  order  in  velocity  (Palczewski  et  ah,  1997;  Palczewski  and 
Schneider,  1998;  Panferov  and  Heintz,  2002). 

Another  important  class  of  numerical  methods  is  based  on  the  use  of  spectral  tech¬ 
niques  in  the  velocity  space.  The  methods  were  first  derived  in  Pareschi  and  Perthame 
(1996),  inspired  by  previous  works  on  the  use  of  Fourier  transform  techniques  (see  Bobylev 
(1988)  for  instance).  The  numerical  method  is  based  on  approximating  the  distribution 
function  by  a  periodic  function  in  velocity  space,  and  on  its  representation  by  Fourier 
series.  The  resulting  scheme  can  be  evaluated  with  a  computational  cost  of  0(n2). 

The  method  was  further  developed  in  Pareschi  and  Russo  (2000b, c)  where  evolu¬ 
tion  equations  for  the  Fourier  modes  were  explicitly  derived  and  spectral  accuracy  of  the 
method  was  proved.  Strictly  speaking  these  methods  are  not  conservative,  since  they 
preserve  mass,  whereas  momentum  and  energy  are  approximated  with  spectral  accuracy. 
This  trade  off  between  accuracy  and  conservations  seems  to  be  an  unavoidable  compromise 
in  the  development  of  numerical  schemes  for  the  Boltzmann  equation. 

Recently  in  Mouhot  and  Pareschi  (2006,  2004);  Filbet  et  al.  (2006),  using  a  suitable 
representation  of  the  collision  operator,  the  computational  cost  of  spectral  methods  has 
been  reduced  from  0(n2)  to  0(n  log2  n)  without  loosing  the  spectral  accuracy  thus  making 
the  methods  competitive  with  Monte  Carlo.  These  fast  algorithms  are  restricted  to  a 
certain  class  of  particle  interactions  including  pseudo-Maxwell  molecules  (for  dv  =  2)  and 
hard  spheres  (for  dv  =  3).  This  kind  of  approach  has  been  extended  recently  to  construct 
fast  algorithms  for  DVM  models  (Mouhot  and  Pareschi,  2011).  Another  class  of  fast 
solvers  for  the  case  of  radially  symmetric  distribution  functions  has  been  constructed  in 
Markowich  and  Pareschi  (2005). 

We  recall  here  that  the  spectral  method  has  been  applied  also  to  non  homogeneous 
situations  (Filbet  and  Russo,  2003),  to  the  Landau  equation  (Filbet  and  Pareschi,  2003; 
Pareschi  et  al.,  2000,  2003),  where  fast  algorithms  can  be  readily  derived,  to  the  case  of 
granular  gases  (Naldi  et  ah,  2003;  Filbet  et  ah,  2005)  and  more  recently  to  the  case  of 
a  quantum  gas  (Filbet  et  al.,  2011).  Let  us  mention  that  algorithms  based  on  a  Fourier 
transform  approximation  of  the  distribution  function  have  been  constructed  in  Bobylev 
and  Rjasanow  (1997,  1999)  and  more  recently  in  Gamba  and  Tharkabhushanam  (2009, 
2010).  Other  fast  algorithms  for  kinetic  equations  can  be  found  in  Buet  et  al.  (1997); 
Lemon  (1998). 

An  additional  problem  in  the  numerical  solution  of  the  Boltzmann  equation  is  the  time 
step  restriction  in  regions  close  to  the  fluid  dynamic  limit,  i.e.  for  very  small  Knudsen 
number.  In  such  a  cases,  in  fact,  the  mean  collision  time  is  so  small  that  an  explicit  solver 
would  require  a  very  small  time  step,  thus  degrading  the  performance  of  the  method. 

If  the  Knudsen  number  is  uniformly  small,  then  one  usually  does  not  need  a  kinetic 
treatment,  and  the  gas  can  be  very  well  described  by  the  Euler  or  Navier-Stokes  equations. 
There  are  cases,  however,  in  which  the  local  Knudsen  number  varies  over  several  orders 
of  magnitude.  Several  authors  have  tackled  the  problem  in  the  past,  and  there  is  a  large 
literature  on  the  subject  (see  Bennoune  et  al.  (2008);  Caflisch  et  al.  (1997);  Degond  et  al. 
(2005);  Gabetta  et  al.  (1997);  Filbet  and  Jin  (2010);  Tiwari  and  Klar  (1998)  and  the 
references  therein). 

One  possibility  is  to  resort  to  the  so  called  domain  decomposition  techniques:  one 
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could  divide  the  computational  domain  into  two  complementary  domains,  let  us  denote 
them  by  A [ero dynamics]  and  B[oltzmann];  in  A  the  gas  is  well  described  by  either  Euler 
or  Navier-Stokes  equations,  while  in  B  the  gas  needs  a  kinetic  description.  In  some  cases 
the  subdivision  into  such  two  subdomains  may  be  known  a  priori ,  for  example  from  the 
geometry  of  the  domain  of  from  previous  approximate  calculations  of  the  flow,  but  in  most 
cases  the  two  regions  are  themselves  unknown,  and  therefore  they  have  to  be  computed  and 
evolved  as  part  of  the  solution.  Examples  of  Euler-Boltzmann  coupling  can  be  found,  for 
example,  in  Bourgat  et  al.  (1992);  Tiwari  and  Klar  (1998);  Tiwari  (1998),  while  Navier- 
Stokes-Boltzmann  coupling  is  considered  in  Bourgat  et  al.  (1996).  More  recent  works, 
including  hybrid  methods,  can  be  found  in  Schwartzentruber  et  al.  (2007);  Degond  et  al. 
(2007);  Dimarco  and  Pareschi  (2010b). 

The  coupling  between  aerodynamics  and  kinetic  description  is  very  appealing,  because 
one  could  gain  maximum  efficiency  by  treating  with  continuum  equations  the  region  in 
which  a  kinetic  description  is  not  strictly  necessary.  However,  it  presents  several  difficul¬ 
ties,  and  it  is  still  an  open  problem  to  asses  what  is  the  best  coupling  strategy. 

In  these  notes  we  shall  consider  system  which  is  treated  by  kinetic  equations  in  the 
whole  computational  domain.  We  shall  discuss  what  strategies  can  be  used  to  treat 
regions  with  small  Knudsen  number  avoiding  unnecessary  restrictions  on  the  time  step. 
When  the  mean  collision  time  is  much  smaller  that  typical  length  scales  of  the  problem, 
the  Boltzmann  equation  becomes  stiff.  Usually  stiff  systems  governed  by  time  dependent 
equations  can  be  effectively  treated  by  implicit  schemes.  The  interested  reader  may  consult 
the  monograph  by  Hairer  and  Wanner  on  numerical  solution  of  stiff  problems  (Hairer  and 
Wanner,  1996).  A  direct  time  discretization  of  the  Boltzmann  equation  seems  not  possible 
in  such  stiff  regimes  due  to  the  high  dimensionality  and  the  nonlinearity  of  the  collision 
operator  which  makes  unpractical  the  use  of  implicit  solvers. 

Here  we  present  two  classes  of  methods  which  avoid  the  solution  of  systems  of  non¬ 
linear  equations.  The  Erst  one  is  based  on  operator  splitting  approach  combined  with 
the  use  of  exponential  techniques,  and  will  be  applied  to  the  Boltzmann  equation  itself 
(see  Gabetta  et  al.  (1997);  Dimarco  and  Pareschi  (2010a)).  The  second  class  is  based  on 
non  splitting  approaches.  We  focus  our  attention  on  implicit  semilagrangian  schemes  (see 
Santagati  (2007);  Russo  and  Santagati  (2011))  and  Implicit-Explicit  Runge-Kutta  meth¬ 
ods  (see  Pareschi  and  Russo  (2005);  Pieraccini  and  Puppo  (2007));  such  methods  will 
be  constructed  for  the  BGK  model  of  the  Boltzmann  equation  and  then  their  extension 
to  the  full  Boltzmann  equation  discussed  (Filbet  and  Jin,  2010;  Dimarco  and  Pareschi, 
2011).  We  refer  the  reader  to  Bennoune  et  al.  (2008)  for  a  related  approach. 

In  addition  to  other  deterministic  methods  that  will  not  be  discussed  in  the  notes,  such 
as  finite  difference  methods  (Ohwada,  1993;  Sone  et  ah,  1989),  there  are  several  important 
aspects  that  we  shall  not  be  able  to  discuss  in  the  present  paper,  namely: 

•  Numerical  treatment  of  boundary  conditions.  In  most  applications  reflective  or 
absorbing  boundary  conditions  should  be  considered,  or  a  suitable  combination  of 
the  two.  However  a  detailed  treatment  of  boundary  conditions  depends  on  the 
geometry  of  the  domain,  and  on  the  detail  of  the  space  discretization  (Carrillo 
et  ah,  2006). 

•  Large  departure  from  (global)  equilibrium.  One  of  the  weak  points  of  deterministic 
methods  based  on  a  fixed  grid  in  velocity  space  is  that  a  huge  grid  in  velocity  would 
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be  necessary  to  treat  flows  with  a  large  variation  of  macroscopic  mean  velocity  and 
temperature.  Effective  schemes  for  the  treatment  of  such  cases  have  to  be  based  on 
the  use  of  grid  in  velocity  domain  which  change  from  one  space  location  to  another 
(Filbet  and  Russo,  2006;  Heintz  et  al.,  2008). 

•  Multiple  space  regimes.  Most  realistic  flows  present  strong  non  homogeneities,  with 
small  regions  of  large  gradients  in  the  moments.  Even  the  region  requiring  kinetic 
treatment  may  present  several  space  scales.  An  effective  treatment  of  such  problems 
would  require  adaptivity  of  the  grid  in  space  (Cai  and  Li,  2010). 

•  Effective  techniques  for  stationary  flows.  These  notes  provide  a  review  of  deter¬ 
ministic  methods  for  the  time  dependent  Boltzmann  equation.  In  fact,  thanks  to 
averaging  procedures,  the  stationary  case  is  usually  computed  efficiently  with  Monte 
Carlo  methods  (Bird,  1994).  However,  one  may  be  interested  in  accurate  determin¬ 
istic  computations  of  the  stationary  solutions,  which  may  be  treated  by  schemes 
aimed  to  capture  the  stationary  state  (Greenberg  and  Leroux,  1996;  Botchorishvili 
et  al.,  2003). 

•  Diffusion  limits.  In  several  circumstances  the  space-time  scaling  of  the  kinetic  equa¬ 
tion  leads  asymptotically  to  the  corresponding  diffusion  system.  These  problems 
present  additional  difficulties  compare  to  the  standard  fluid  scaling  since  also  the 
transport  terms  are  strongly  stiff  (Jin  et  al.,  2000;  Lernou  and  Mieussens,  2008). 


1  The  Boltzmann  equation 

1.1  The  model 

The  model  is  characterized  by  a  density  function  f(x,v,t)  describing  the  time  evolution 
of  a  monoatomic  rarefied  gas  of  particles  which  move  with  velocity  v  e  JR 3  in  the  position 
x  G  H  C  JR 3  at  time  t  >  0  which  satisfies  the  Boltzmann  equation  (Cercignani,  1988; 
Cercignani  et  al.,  1994) 

%  +  v-VJ=-Q(f,f),  (1) 

at  e 

with  initial  data 

f{x,v,0)  =  f0(x,v).  (2) 

The  parameter  £  >  0  is  called  Knudsen  number  and  is  proportional  to  the  mean  free 
path  between  collisions.  The  bilinear  collision  operator  Q(f,  f)  which  describes  the  binary 
collisions  of  the  particles  acts  over  the  velocity  variable  only 

Q(f,f)(v)=  f  [  B(v,v*,u)[f(v')f(v*)- f(v)f(v*)]dudv*.  (3) 

J  m3  J  s2 

In  the  above  expression,  u  is  a  unit  vector  of  the  sphere  §2  and  (vf,  u*)  represent  the 
collisional  velocities  associated  with  (v,v*).  The  collisional  velocities  satisfy  microscopic 
momentum  and  energy  conservation 

v'  +  v*  —  v  +  v*,  \v'\2  +  |t/  j2  =  H2  +  |u*|2.  (4) 
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The  above  system  of  algebraic  equations  has  the  following  parametrized  solution 
1  1 

v  —  g +  \v  ~  v*\u),  v '  =  ~{v  +  u*  —  \v  —  v*\ u)  (5) 

where  v  —  u*  is  the  relative  velocity. 

The  collision  kernel  B(v,  v*,uu)  is  a  nonnegative  function  which  characterizes  the  details 
of  the  binary  interactions  and  depends  only  on  |u  —  v*|  and  the  scattering  angle  9  between 
relative  velocities  v  —  u*  and  v'  —  v*  =  \v  —  u*| u 

a  _  {v  -  v*)  '  ^ 

V  -  u* 


The  kernel  has  the  from 


B(v,v*,u>)  =  \v  -  u*|cr(|u  -  u*|,cos0),  (6) 

where  the  function  a  is  the  scattering  cross-section. 


Example  1 

•  In  the  hard  sphere  model  the  particles  are  assumed  to  be  ideally  elastic  spheres  of 
diameter  d  >  0  and  thus 


<t(|u  —  u*|,  cos  9)  —  — ,  B(v,v*,u)  =  —  \v  —  v*\,  (7) 

since  the  total  cross  section  is  7r d2  =  47r(d2/4). 

•  In  the  case  of  inverse  k-th  power  forces  between  particles,  the  kernel  has  the  form 

a(\v  —  u*|,  cos  9)  =  ba(cos9)\v  —  B(v,  v*,uf)  =  ba( cos  9)\v  —  u*|a,  (8) 

with  a  =  (k  —  5)/(k  —  1).  For  k  >  5  we  have  hard  potentials,  for  k  <  5  we  have 
soft  potentials. 

•  The  special  situation  k  =  5  gives  the  Maxellian  model  with 

B(v,v*,u)  —  b0(cos9).  (9) 


•  For  numerical  purposes,  a  widely  used  model  is  the  Variable  Hard  Sphere  (VHS) 
model,  corresponding  to  ba(cos9)  =  Ca,  where  Ca  is  a  positive  constant,  and  hence 

a(\v  -  d*|,cos0)  =  Ca\v  -  u*|a_1,  B(v,v*,u)  =  Ca \v-v*\a.  (10) 
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1.1  The  model 


Figure  1:  The  collision  sphere 


The  collision  integral  Q(f,  f )  can  be  written  in  different  equivalent  forms,  according 
to  the  parametrization  used  for  the  collisional  velocities.  Using  the  identity 


Jju  ■  n)+4>{n{u  •  n))  dn  =  H  j  cf>  I  j  du 

obtained  by  the  transformation  uj  —  e  —  2(e  •  n)n,  we  get  the  frequently  used  form 


(11) 


with 


and 


Q(f,f)(v)  =  B(v,v*,u)[f(v')f(v'J  -  f(v)f(v*)}dudv* 

J m3  Js2 


v'  =  v  —  ((v  —  iu)  ■  u)u,  w'  =  u*  +  ((v  —  u*)  •  u)oj, 


B(v,v*,uj)  =  2\v  —  d*||  cosd|cr(|n  —  u*|,  1  —  2|  cos 6*1). 
The  hard  sphere  case  corresponds  to 

d2 

B(v,v*,uj)  =  —\v  -  u*||  cos0|, 
whereas  the  Maxwellian  molecules  case  gives 

B(v,v*,uj)  =  2|  cos#|f>0(cos#). 


(12) 

(13) 

(14) 

(15) 

(16) 


Remark  1 

For  the  Maxwellian  case  the  collision  kernel  B(v,v*,lj )  is  independent  of  the  rela¬ 
tive  velocity.  This  case  has  been  widely  studied  theoretically,  in  particular  exact  analytic 
solutions  can  be  found  in  the  space  homogeneous  case  where  f  =  f(v,t )  (Bobylev,  1975). 
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A  simplified  one- dimensional  space  homogeneous  Maxwell  model  is  given  by  the  Kac 
equation  (Kac,  1957).  It  reads 

%  =  JMfQ  ~  f(v)f(v*)}  dedv*  (17) 

where  the  collisional  velocities  are  characterized  by  rotations  in  the  collisional  plane 

v*  =  v  cos  9  —  v*  sin  9 ,  v*  =  v  sin  9  +  v*  cos  9.  (18) 


For  this  model  we  have  only  microscopic  conservation  of  energy  ( v ')2  +  (v()2  =  v2  +  v2. 


1.2  Physical  properties 

During  the  evolution  process,  the  collision  operator  preserves  mass,  momentum  and  en¬ 
ergy,  i.e., 

[  Q(fJ)</>(v)dv  =  0,  <j>(v)  =  l,vx,vv,vz,\v\2,  (19) 

Jn 3 

and  in  addition  it  satisfies  Boltzmann’s  well-known  //-theorem 


[  Q(f,f)\n(f(v))dv  <  0-  (20) 

Jm3 

The  above  properties  are  a  consequence  of  the  following  identity  that  can  be  easily  proved 
for  any  test  function  0(u) 


dv 


T  f  [  B(v,  v*,u)  [/'/*  -  //*]  W  +  0*  -  0  -  0*]  du  dry  dv. 
4  J]R6  J s2 


where  we  have  omitted  the  explicit  dependence  from  v,  u*,  v',  v*  to  simplify  the  expression. 

In  order  to  prove  this  identity  we  used  the  micro- reversibility  property  B{y,v*,u)  = 
B(v*,v,u )  and  the  fact  that  the  Jacobian  of  the  transformation  (v,v*)  fG  (v',v*)  is  equal 
to  1. 

A  function  0  such  that 


fiW)  +  0K)  ~  &{v)  ~  Hv*)  =  0 

is  called  a  collision  invariant.  It  can  be  shown  that  a  continuous  function  0  is  a  collision 
invariant  if  and  only  if  0  e  span{l,u,  |u|2}  or  equivalentely 

0(u)  —  a  +  b  ■  v  +  c\v\2,  a,c  E  IR,  be  IR3. 

Assuming  /  strictly  positive,  for  0(u)  =  In (f(v))  we  obtain 

[  Q(f,f)\n(f)dv 
JlR3 

=  -  J  [  [  B(v,v*,u)[ff*  -  //*][ln(//)  +  In (/')  -  In (/)  -  ln(/*)]  du  dv*  dv 

4  J  w  J  s2 

=  f  [  B(v,v*,u)[ffl  -  //*]  In  dudv*dv  <  0, 

4  0JR6  J S2  \JJ*J 
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1.3  Fluid  limit 


since  the  function  z(x,  y)  —  (x  —  y )  In (x/y)  >  0  and  z(x,  y)  —  0  only  if  x  =  y. 

In  particular,  the  equality  holds  only  if  ln(/)  is  a  collision  invariant,  that  is 

/  =  exp(a  +  b  ■  v  +  c|u|2),  c  <  0. 

If  we  define  the  density,  mean  velocity  and  temperature  of  the  gas  by 

P=  [  fdv,  u  =  -  [  vfdv,  T  =  ^—  j  [ v-uffdv ,  (21) 

Jr3  P  J r3  6KP  Jr3 

we  obtain  that  the  distribution  function  has  the  form  of  a  locally  Maxwellian  distribution 

/M)  =  M(p,u,T)(v,t)  =  {2JT)3/2  exp  (-U^L)  • 

The  constant  R  =  Kb/tu  is  called  the  gas  constant,  Kb  is  the  Boltzmann  constant  and  m 
the  mass  of  a  particle.  Boltzmann’s  //-theorem  implies  that  any  equilibrium  distribution 
function,  i.e.  any  function  /  for  which  Q(f,  f)  =  0,  has  the  form  of  a  locally  Maxwellian 
distribution. 

If  we  define  the  H-function 

//(/)=/  f\n(f)dv, 

Jr3 

we  obtain  immediately  the  inequality 

«(/,/)  ln(/)  dv  <  0.  (22) 

at  jM  3 

Thus  the  //-function  is  monotonically  decreasing  until  /  reaches  the  equilibrium  Maxwellian 
state  for  which  we  have 

H(M)  =  P  (hl  ((2t tRT)3/2)  ~  2)  ' 

1.3  Fluid  limit 

If  we  multiply  the  Boltzmann  equation  by  its  collision  invariants  and  integrate  the  result 
in  velocity  space  we  obtain 

4  dv  +  (  [  vfHv)  dv 

9t  Jr3  \Jr3 

These  equations  describe  the  balance  of  mass,  momentum  and  energy.  The  system  of  five 
equations  is  not  closed  since  it  involves  higher  order  moments  of  the  distribution  function 
/• 

As  £  — »  0,  from  (1)  we  have  formally  Q(f,  f)  — »  0,  and  thus  /  approaches  the  local 
Maxwellian.  In  this  case  the  higher  order  moments  of  the  distribution  function  can  be 


^  =  0,  4>{v)  =  l,v1,v2,v3,\v 
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computed  as  function  of  p,  u,  and  T  and  we  obtain  the  closed  system  of  compressible 
Euler  equations 

^  +  Vx  •  (pu)  =  0 
dpu  „ 

+  V;c  •  (pu  <g>  u  +  p)  =  0 
dE 

—  +  Vx  ■  ( Eu  +  pu)  =  0 
P  =  pT,  E  =  | pT  +  ^pw2 

where  p  is  the  gas  pressure  and  <8)  denotes  the  tensor  product. 

The  rigorous  passage  from  the  Boltzmann  equation  to  the  compressible  Euler  equations 
has  been  investigated  by  several  authors.  Among  them  we  mention  references  Caflisch 
(1980);  Nishida  (1978).  Higher  order  fluid  models,  such  as  the  Navier-Stokes  model, 
can  be  considered  using  the  Chapmann-Enskog  and  the  Hilbert  expansions.  We  refer 
to  Levermore  (1996)  for  a  mathematical  setting  of  the  problem  and  to  Golse  and  Saint- 
Raymond  (2004)  for  recent  theoretical  results. 


1.4  Boundary  conditions 

The  Boltzmann  equation  is  complemented  with  the  boundary  conditions  in  space  for 
v  ■  n  >  0  and  x  G  (90,  where  n  denotes  the  unit  normal,  pointing  inside  the  domain  ft 
Usually  the  boundary  represents  the  surface  of  a  solid  object  (an  obstacle  or  a  container). 
The  particles  of  the  gas  that  hit  the  surface  interact  with  the  atoms  of  the  object  and  are 
reflected  back  into  the  domain  U. 

Mathematically,  such  boundary  conditions  are  modelled  by  an  expression  of  the  form 
(Cercignani,  1988) 

u*  -n(x)\K(v*  v,x,t)f(x,v*,t)  dv*.  (23) 

This  is  the  so-called  reflective  boundary  condition  on  dfi. 

The  ingoing  flux  is  defined  in  terms  of  the  outgoing  flux  modified  by  a  given  boundary 
kernel  K .  This  boundary  kernel  is  such  that  positivity  and  mass  conservation  at  the 


\v  ■  n\f(x,  v,  t)  = 


'  v*-n<0 


Figure  2:  Reflection  and  diffusion  at  the  solid  boundary 
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1.4  Boundary  conditions 


boundaries  are  guaranteed, 


K(v*  — »  u,  x,  t)  >  0,  /  A (u*  — *  u,  x,  t )  dv  =  1. 

J  v-n(x)>  0 

Commonly  used  reflecting  boundary  conditions  are  the  so-called  Maxwell’s  conditions. 
From  a  physical  point  of  view,  one  assumes  that  a  fraction  a  of  molecules  is  absorbed  by 
the  wall  and  then  re-emitted  with  the  velocities  corresponding  to  those  in  a  still  gas  at 
the  temperature  of  the  wall,  while  the  remaining  fraction  (1  —  a)  is  specularly  reflected. 

This  is  equivalent  to  impose  for  the  ingoing  velocities 

f(x,v,t)  =  (1  -  a)Rf(x,v,t)  +  aMf(x,v,t),  (24) 

in  which  x  G  <9h2,  v  ■  n{x )  >  0.  The  coefficient  a,  with  0  <  a  <  1,  is  called  the 
accommodation  coefficient  and 

Rf(x,v,t)  =  f(x,v  -*  2n(n  ■  v),t),  Mf(x,  v,  t)  =  n(x,  t)Mw[v).  (25) 


If  we  denote  by  Tw  the  temperature  of  the  solid  boundary,  Mw  is  given  by 


Mw(v)  =  exp(- 


2  RTwn 

and  the  value  of  /j  is  determined  by  mass  conservation  at  the  surface  of  the  wall 


/x(x,t)  /  Mw(y)\v  ■  n\dv  —  /  f(x,v,t)\v-n\dv. 

J  v-n> 0  J  v-n< 0 


(26) 


For  a  =  0  (specular  reflection)  the  re-emitted  molecules  have  the  same  flow  of  mass, 
temperature  and  tangential  momentum  of  the  incoming  molecules,  while  for  a  =  1  (full 
accommodation)  the  re-emitted  molecules  have  completely  lost  memory  of  the  incoming 
molecules,  except  for  conservation  of  the  number  of  molecules. 

More  complex  boundary  conditions  for  rarefied  gas  dynamics  (RGD)  can  be  imposed 
using  the  boundary  conditions  of  Cercignani  and  Larnpis  (Cercignani  and  Lampis,  1971). 
These  can  be  written  as 


f{x,v,t)  =  J  P(v,v')f(x,v',t)dv'  (27) 

where 

P(v,  v  )  =  {2v / a)I' {2{1  —  a)l / 2vv' / a)  exp  ( v 2  —  (1  —  a)vr2)/a  (28) 

in  which  v  and  v'  are  the  normal  components  of  the  outgoing  and  incoming  velocities 
respectively,  and  I'  is  the  modified  Bessel  function.  This  satisfies  the  reciprocity  (detailed 
balance)  condition 

vP(-v',  -v)Mw(v)  =  -v'P(v,  v')Mw{v').  (29) 

A  consequence  of  the  reciprocity  condition  is  that  the  Maxwellian  distribution  Mw  is 
preserved  by  this  boundary  condition. 
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In  the  case  of  inflow  boundary  conditions,  one  assumes  that  the  distribution  function 
of  the  particles  entering  the  domain  is  known,  i.e. 

f(x,v,t)—g(v,t),  x  G  <90,  v  ■  n  >  0, 

A  typical  example  of  such  condition  is  used  in  shock  wave  calculations,  where  one  as¬ 
sumes  that  the  distribution  function  at  the  boundary  of  the  computational  domain  is  a 
Maxwellian  M (v)  and  that  the  incoming  flux  of  particles  at  the  boundary  is  distributed 
according  to  the  Maxwellian  flux  (v  ■  n)M{y),  v  ■  n  >  0. 

1.5  Other  collision  operators 

1.5.1  BGK  models 

A  simplified  model  Boltzmann  equation  is  the  so-called  BGK  model  introduced  by  Bhat- 
nagar  et  al.  (1954).  In  this  model  the  collision  operator  is  replaced  by  a  relaxation  operator 
of  the  form 

QBGK(f,f)(v)  =  -(M{f]~f),  (30) 

T 

where  M[f]  =  M(v,  {p,u,T})  is  the  local  Maxwellian  computed  by  the  moments  of  the 
distribution  function  / 

M(v-  {„,  u,T})=  (2jrJ£,)3/2  exp  (-K=jA)  •  (31) 

where  p  =  p(x,t),  u  =  u(x,t )  and  T  =  T(x,t)  denote  the  macroscopic  fields,  namely: 
density,  mean  velocity  and  temperature,  corresponding  to  the  function  /. 

Conservation  of  mass,  momentum  and  energy  as  well  as  Boltzmann  H-theorem  are 
readily  satisfied.  The  equilibrium  solutions  are  clearly  Maxwellians 

Qbgk(/,/)  =  0^/  =  M[/]. 

The  relaxation  time  r  is  in  general  inversely  proportional  to  the  density,  and  depends  on 
the  temperature: 

r-1  =  A(T)p 

Numerical  computations,  as  well  as  the  analytic  theory,  for  such  model  are  much  simpler 
then  for  the  full  Boltzmann  equation. 

Furthermore,  as  a  consequence  of  conservation  of  mass,  momentum  and  energy,  in  the 
fluid  dynamic  limit  the  moments  (p, 

rhou ,  and  E ),  i.e.  mass  density,  momentum  density,  and  total  energy  density,  satisfy  the 
compressible  Euler  equations  for  a  monoatomic  gas,  therefore  the  model  describes  the 
correct  fluid  dynamic  limit. 

But  in  the  Chapman-Enskog  expansion,  the  transport  coefficients  obtained  at  the 
Navier-Stokes  level  are  not  satisfactory.  The  relaxation  time  could  be  adjusted  so  that  at 
the  Navier-Stokes  level  the  model  provides  the  correct  value  of  one  transport  coefficient, 
say  the  viscosity.  However,  the  Prandtl  number  Pr  (the  ratio  between  heat  conductivity 
and  viscosity)  is  equal  to  1.  For  most  gases,  we  have  Pr  <  1.  In  particular,  the  hard-sphere 
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model  for  a  monoatomic  gas  leads  to  a  Prandtl  number  very  close  to  2/3,  therefore  only 
one  transport  coefficient  can  be  correct,  but  not  both.  The  correct  Prandtl  number  can 
be  recovered  using  more  sophisticated  BGK  models,  as  the  velocity  dependent  collision 
frequency  BGK  models  and  the  Ellipsoidal  Statistical  BGK  (ES-BGK)  models  (Bouchut 
and  Perthame,  1993;  Holway,  1966). 

1.5.2  Landau  models 

The  Landau  model  (Landau,  1981)  is  a  common  kinetic  model  in  plasma  physics  charac¬ 
terized  by  the  following  collision  operator 

Ql(I,  f)(v)  =  Vv  ■  [  A(v-  u*)[/(v*)V„/(u)  -  f(v)VvJ{v*)]  du* 

J  Md 

where  A(z)  =  \I/(|;z|)n(z)  is  a  d  x  d  nonnegative  symmetric  matrix  and  n(^)  =  (ir ’ij(z)) 
is  the  orthogonal  projection  upon  the  space  orthogonal  to  z, 

-  |j)  • 

We  have  H/(|^|)  =  A|z|"+2  for  inverse-power  laws,  with  a  >  —3  and  A  >  0.  The  case 
oi  =  —3  is  the  so-called  Coulombian  case,  of  primary  importance  for  applications.  In 
such  case  the  Boltzmann  collision  operator  has  no  meaning,  due  to  the  divergence  of  the 
integral,  even  for  smooth  functions  (a  cut-off  angular  approximation  is  then  used  and  the 
Landau  equation  can  be  derived  in  the  so  called  grazing  collision  limit  (Villani,  2002)). 

Since  conservation  of  mass,  momentum,  and  energy,  as  well  as  H-theorem  for  the 
entropy  are  satisfied,  equilibrium  states  are  Maxwellians. 

1.5.3  Additional  models 

•  Enskog  model:  takes  into  account  the  nonlocality  of  the  interactions  induced  by 
the  diameter  of  the  interacting  spheres  (accurately  describes  the  behavior  of  dense 
gases).  The  collision  operator  is  delocalized  in  space  (regularization  effect). 

•  Quantum-Boltzmann  models:  the  nonlinear  interactions  /'/*  —  //*  is  replaced  by 

/7'(1  ±  /)(1  ±  /*)  -  //*(1  ±  /')(1  ±  /*)• 

The  minus  sign  corresponds  corresponds  to  fermions  (such  as  alectrons),  and  the 
plus  sign  to  bosons  (such  as  photons).  The  collision  operator  are  called  Pauli  oper¬ 
ator  and  Bose-Einstein  operator  respectively. 

•  Semiconductor- Boltzmann  models:  the  linear  Boltzmann  equation  for  semiconduc¬ 
tor  devices  has  the  form 

Qs(f,  M)  —  j  <r(v,v*){M(v)f{v*)  -  M(v*)f(v)}dv*, 

where  M  is  the  normalized  equilibrium  distribution  (Maxwellian,  Fermi-Dirac)  at 
the  temperature  6  of  the  lattice.  The  function  a(v,v*)  describes  the  interaction  of 
carriers  with  phonons. 
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•  Granular  gas  models:  particles  undergo  inelastic  collisions.  Energy  is  dissipated 
by  the  model  and  the  steady  states  are  Dirac  delta  function  centered  in  the  mean 
velocity. 

More  recently  kinetic  modelling  has  been  applied  to  new  fields  as  vehicular  traffic  flows, 
biomathematics  (chemotaxis,  inhalation  of  sprays,  flocking),  finance  (modelling  income 
distributions,  price  formation),  coagulation- fragmentation  processes,  supply  chains,  and 
so  on. 

1.6  The  splitting  approach 

The  most  common  approach  to  solve  numerically  the  full  Boltzmann  equation  is  based  on 
an  operator  splitting  (Desvillettes  and  Mischlcr,  1996). 

The  solution  in  one  time  step  At  may  be  obtained  by  the  sequence  of  two  steps.  First 
integrate  the  space  homogeneous  equation  for  all 

%  =  ?W,/),  (32) 

f(x,v,  0)  =  fo(x,v), 

for  a  time  step  At  ( collision  step )  to  obtain  /  =  Ca*(/o),  and  then  the  transport  equation 
using  the  output  of  the  previous  step  as  initial  condition, 

^+v-Vxf  =  0,  (33) 

f(x,v,  0)  =  f(x,v,At). 

for  a  time  step  At  ( transport  step )  to  get  /  =  7a*(/)  =  7At(CAt(/o))- 

After  computing  an  approximation  of  the  solution  at  time  At,  the  process  may  be 
iterated  to  obtain  the  numerical  solution  at  later  times.  Although  this  splitting  scheme 
(simple  splitting)  described  above  is  first  order  accurate  in  time  it  is  very  popular  because 
it  has  several  nice  properties. 

•  The  collision  step  acts  only  on  v  whereas  the  transport  step  acts  on  x.  This  makes 
the  implementation  of  the  resulting  scheme  simpler  (it  allows  the  use  of  any  existing 
code  designed  to  solve  the  free  transport  equation)  and  highly  parallelizable. 

•  It  makes  simpler  to  design  schemes  which  preserves  the  physical  properties  of  the 
equation  (conservations,  positivity,  H-theorem),  since  these  properties  essentially 
depends  on  the  treatment  of  the  collision  step. 

It  is  then  clear,  that  after  this  splitting  almost  all  the  main  numerical  difficulties  are  con¬ 
tained  in  the  collision  step.  The  discretization  of  the  resulting  equations  can  be  performed 
in  a  variety  of  ways  (finite  volume,  finite  difference,  Monte  Carlo  methods  and  so  on).  The 
choice  of  the  discretization  mainly  depends  on  the  method  that  is  used  for  the  solution 
of  the  space  homogeneous  Boltzmann  equation.  Higher  order  splitting  formulas  can  be 
derived  in  different  ways  (see  Hairer  et  al.  (2002)).  For  example  the  well-known  second 
order  Strang  splitting  (Strang,  1968)  can  be  written  as 

At/2  (Taj  (C At/2  (/o)  )  )  ■  (34) 
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1.7  Asymptotic  preserving  methods 


Unfortunately  for  splitting  methods  of  order  higher  then  two  it  can  be  shown  that  it’s 
impossible  to  avoid  negative  time  steps  both  in  the  transport  as  well  as  in  the  collision 
(Hairer  et  ah,  2002).  As  an  example  we  report  here  a  symmetric  fourth  order  formula 
(McLachlan,  1995) 

'Ta1At{Cb1At{'Ta2At{Cb2At{'Ta3At{Cb2At{Ta2At(Cb1At(Ta1At{fo))))))))),  (35) 

where 

61  =  ^-,  b2  =  i  -  61  «  -0.045,  (36) 

.  =  ^,0.169,  a2  =  tr^(12- V471)  ~  -0.299,  a3  =  l-2(a1+a2).  (37) 

Higher  order  formulas  which  avoid  negative  time  stepping  can  be  obtained  as  suitable 
combination  of  splitting  steps  (Dia  and  Schatzman,  1996).  For  example  a  third  order 
approximation  is  given  by 

2  [7At/2(CAt(7At/2(/o)))  +  CAt/2(TAt(CAt/2(fo)))]  ~  g  [7At(dAt  (/o))  +  (7a*(/o))]  j  (38) 

which  corresponds  to  take  a  combination  of  symmetrized  Strang  and  first  order  splitting, 
whereas  a  fourth  order  scheme  reads 

-CAt/4(TAt/2(CAt/2{TAt./2(CAt/4:(fo)))))  ~  “C At/2(7At(^At/2(fo )))•  (39) 

Clearly  all  the  above  splitting  methods  admit  the  symmetric  formulation  obtained  by 
switching  the  transport  and  the  collision  operators.  Note,  however,  that  the  appearance 
of  negative  coefficients  or  negative  time  steps  in  high  order  formulas  may  lead  to  some 
drawbacks  in  practical  applications  like  the  lack  of  positivity  of  the  solution  which  makes 
very  difficult  their  use  in  Monte  Carlo  schemes. 

1.7  Asymptotic  preserving  methods 

Even  if  it  is  difficult  to  give  a  rigorous  definition  of  asymptotic  preserving  scheme  since 
the  concept  has  been  used  for  a  long  time  in  the  physics  and  mathematics  literature  and 
may  refer  to  different  discretization  parameters  (see  Figure  3),  here  following  Jin  (1995) 
and  Pareschi  and  Russo  (2005),  we  formalize  this  notion  for  the  time  discretization  of 
equation  (1). 

Definition  1  A  consistent  time  discretization  method  for  (1)  of  stepsize  At  is  asymptotic 
preserving  (AP)  if,  independently  of  the  initial  data  (2)  and  of  the  stepsize  At,  in  the  limit 
£  — V  0  becomes  a  consistent  time  discretization  method  for  the  reduced  system  (23). 

Note  that  this  definition  does  not  imply  that  the  scheme  preserves  the  order  of  accuracy 
in  t  in  the  stiff  limit  £  — >  0. 

In  the  case  of  operator  splitting  we  can  reformulate  the  asymptotic  preserving  property 
and  prove  that 
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£  — y  0 

pe _ ^  p  0 


h  ->•  0 


0 


£  — >•  0 


Figure  3:  The  AP  diagram.  Here  P£  is  the  original  singular  perturbation  problem  and 
Pf  its  numerical  approximation  characterized  by  a  discretization  parameter  h.  The  AP 
property  corresponds  to  the  request  that  P f  is  consistent  with  Ps  as  £  — >  0  independently 
of  h. 

Proposition  1  A  sufficient  condition  for  a  consistent  time  discretization  method  of  step- 
size  At  applied  to  the  operator  splitting  approximation  of  (1),  given  by  (32)- (33),  to  be 
AP  is  that  the  time  discretization  of  step  (32),  independently  of  the  initial  data  (2)  and 
of  the  stepsize  At,  in  the  limit  e  — *  0  projects  the  solution  f  over  the  local  Maxwellian 
equilibrium  M (p0,  u0,  T0 ) . 

The  proof  of  the  above  proposition  is  an  immediate  consequence  of  the  fact  that  as 
£  — »  0  step  (32)  degenerates  into  the  projection  CAtffo )  =  M(p0,uo,T0 )  which  coupled 
with  the  transport  step  (33)  originates  a  so-called  kinetic  approximation  (Coron  and 
Perthame,  1991)  to  the  Euler  equation  (23)  given  by  TAt(M (p0,  uq,T0)) .  We  omit  further 
details. 

In  other  words,  Proposition  1  states  that  if  the  relaxation  step  (32)  is  AP  then  the 
whole  splitting  (32)-(33)  is  AP.  Analogous  results  hold  true  for  the  higher  order  splitting 
methods  (34),  (35),  (38)  and  (39).  Let  us  point  out  that  degradation  to  first  order 
accuracy  when  £  — )•  0  is  observed  for  most  splitting  methods  like  the  one  reviewed  here.  A 
possible  way  to  overcome  this  drawback  is  based  on  the  use  of  Implicit-Explicit  (IMEX) 
Runge-Kutta  methods  for  the  full  problem  (1).  We  will  discuss  this  aspect  in  Section  4.3. 

For  the  sake  of  completeness  we  finally  introduce  the  notion  of  entropic  stability, 
namely  schemes  that  preserve  at  a  discrete  level  the  entropy  inequality  (22).  Let  us 
denote  by  fn,  n  >  1  the  numerical  solution  at  t  =  nAt  obtained  with  a  given  time 
discretization  method  applied  to  (32)  with  initial  data  /0. 

Definition  2  A  time  discretization  method  for  (32)  is  called  unconditionally  entropic  if 
H(fn+1 )  <  H(fn),  where  H(f)  =  fR3  flogfdv,  independently  of  the  step  size  At. 

As  pointed  out  in  Dimarco  and  Pareschi  (2010a)  and  Higueras  (2005),  except  for  first 
order  implicit  Euler,  the  entropy  inequality  is  not  satisfied  by  high  order  implicit  Runge- 
Kutta  schemes  applied  to  (32)  unless  a  suitable  time  step  restriction  is  considered.  At 
variance  exponential  methods  permits  to  construct  unconditionally  entropic  methods  at 
any  order  of  accuracy  (Gabetta  et  ah,  1997). 
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In  the  next  section  we  will  focus  on  the  solution  to  the  space  homogeneous  Boltzmann 
equation  (32).  It  is  clear,  in  fact,  that  most  computational  challenges  related  to  the 
behavior  of  the  full  equation  depend  on  the  way  we  approximate  the  collision  operator. 


2  Fast  Boltzmann  solvers 

In  this  section  we  shall  approximate  the  collision  operator  starting  from  a  representation 
which  somehow  conserves  more  symmetries  of  the  collision  operator  when  one  truncates 
it  in  a  bounded  domain.  This  representation  was  used  in  Bobylev  and  Rjasanow  (1997), 
Bobylev  and  Rjasanow  (1999),  Bobylev  and  Rjasanow  (2000),  Ibragimov  and  Rjasanow 
(2002)  and  it’s  close  to  the  classical  Carleman  representation  (cf.  Carleman  (1932)).  As  we 
will  see  it  is  an  essential  step  for  the  derivation  of  fast  algorithms.  The  presentation  here 
follows  the  line  developed  in  Mouhot  and  Pareschi  (2006),  Mouhot  and  Pareschi  (2004), 
Filbet  et  al.  (2006),  Mouhot  and  Pareschi  (2011). 


2.1  Restriction  to  bounded  domains 


The  basic  identity  we  shall  need  is 


1 

2 


F(\u\u  —  u )  du 


_  /  5(2  x  ■  u  +  hi2)  Fix)  dx, 
W~2  V 


(40) 


and  can  be  verified  easily  by  completing  the  square  in  the  delta  Dirac  function,  taking 
the  spherical  coordinate  x  =  ru  and  performing  the  change  of  variable  r2  =  s. 

Setting  u  =  v  —  u*  we  can  write  the  collision  operator  in  the  form 


Q(f,f)(v)=  /  \  I  B(cos6,\u\) 

Jv I  Jue S'*-1 


f(v*  -  n)  f(y  +  -  /(«*)  f(v) 


du  >  dv* 


and  thus  equation  (40)  yields 


Q(f,f)(v)  =  2 


B 


x  ■  u 
\x\\u\ 


u 


U 


—  5{2x  ■  u  +  |x|2) 


ld-2 


f(v*  -  x/2)  f(v  +  x/2)  -  /(u*)  f(v) 


dx  >  dv* 


Now  let  us  make  the  change  of  variable  x  — >•  x/2  in  x  to  get 


«(/,/)( v)=2d+1 


B 


x  ■  u 


XU 


u 


u 


d—  5(4:X  ■  u  +  4|x|2) 


[f(v*  -  x)  f(v  +  x)  -  f(v*)  f(v)}  dxdv* 
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and  then  setting  y  —  v*  —  v  —  x  in  v*  we  obtain 


W, /)(«)  =  2° 


'y£RdJx£Rd  \Fllw 


^5{-4x-y) 


[f(v  +  y)  f(v  +x)-  f(v  +  x  +  y)  f(v)]  dx  dy 
where  now  u  —  —(x  +  y).  Thus  in  the  end  we  have 

W.  f)M  =  2-  £rj  j^d  B  (-^£±d.  |»  +  „|)  5(x  ■  y) 

[f(v  +  y )  f(v  +  x)-  f(v  +  x  +  y)  f(v)}  dx  dy. 

Figure  4  sums  up  the  different  geometrical  quantities  of  the  usual  representation  and  the 
one  we  derived  from  Carleman’s  one. 


\x  +  y\ 


WTW^S(X'V) 


Figure  4:  Geometry  of  the  collision  (v,v*)  O  (v\  u'). 

Now  let  us  consider  the  bounded  domain  VT  =  [— T,  T]d  (0  <  T  <  +oo).  There  are 
two  possibilities  of  truncation  to  reduce  the  collision  process  in  a  box.  From  now  on  let 
us  write 

B(x,y)  =  2 d~lB  (,~~rrr~~^Ti  k  +  J/l)  \x  +  y\~(d~2). 

\  \x\\x  +  y\  ) 

One  can  easily  see  that  on  the  manifold  defined  by  x  ■  y  =  0,  a  simpler  formula  is 

B(x,  y)  =  B(\x\,  b|)  =  2 B  (  == ,  +  \yA  (\x?  +  \y\>)-^ .  (41) 

First  one  can  remove  some  physical  collisions  connecting  with  some  points  out  of  the  box. 
This  is  the  natural  preliminary  stage  for  deriving  conservative  schemes.  In  this  case  there 
is  no  need  for  a  truncation  on  the  modulus  of  x  and  y  since  we  impose  them  to  stay  in 
the  box.  It  yields 


Qtr(f,  f)(v)  = 


B(x,  y)  S(x  ■  y) 

|  v-\-x,  v+y,  v-\-x+y  G  Pj  } 

[f(v  +  y )  f(v  +  x)  -  f(v  +  x  +  y)  f(v)]  dx  dy 
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2.1  Restriction  to  bounded  domains 


defined  for  v  E  T>t ■  One  can  easily  check  that  the  following  weak  form  is  satisfied  by  this 
operator 

[  QtT(f,f)<p(v)dv  =  j  f  f  f  B(x,  y)  S(x  ■  y) 

J  ^  J  J  J  |  v,  x,  y  G  |  v,  v+x,  v+y,  v+x-\-y  G  T>t  | 

f(v  +  x  +  y)  f(v)  [p{v  +  y)  +  <p(v  +  x)  —  (p(v  +  x  +  y)  —  <p{v)\  dv  dx  dy  (42) 

and  this  implies  conservation  of  mass,  momentum  and  energy  as  well  as  the  //-theorem 
on  the  entropy.  Note  that  at  this  level  this  formulation  gives  no  advantage  with  respect 
to  the  usual  one  obtained  from  (3)  by  restricting  v,v*,v',vl  G  T>t-  The  problem  of  this 
truncation  is  that  it  corresponds  to  change  the  collision  kernel  by  adding  some  artificial 
dependence  on  v,v*,v' ,v*.  In  this  way  convolution- like  properties  are  broken. 

A  different  approach  consists  in  periodizing  the  function  /  on  the  domain  T>t ■  Here 
we  have  to  truncate  the  integration  in  x  and  y  since  periodization  would  yield  infinite 

result  if  not.  Thus  we  set  them  to  vary  in  13 r,  the  ball  of  center  0  and  radius  R.  Then  a 

geometrical  argument  (see  Pareschi  and  Russo  (2000b))  shows  that  using  the  periodicity 
of  the  function  it  is  enough  to  take  T  >  (3  +  V2)R/2  to  prevent  intersections  of  the  regions 
where  /  is  different  from  zero. 

The  operator  now  reads 


QR(f,f){v)  =  [  I  B(x,y)S(x-y) 

Jx&Br  Jy&BR 

[f(v  +  y)f(v  +  x)-  f(v  +  X  +  y)f{v)]  dx  dy  (43) 

for  v  e  T>t  (the  expression  for  v  E  is  deduced  by  periodization).  The  interest  of  this 
representation  is  to  preserve  the  real  collision  kernel  and  its  properties. 

By  making  some  translation  changes  of  variable  on  v  (by  x,  y  and  x  +  y),  using  the 
changes  x  — »  —  x  and  y  — >■  —y  and  the  fact  that 

B{-x,  y)  8{-x  ■  y)  =  B(x ,  y)  8(x  ■  y)  =  B(x,  -y)S(x  ■  -y) 

one  can  easily  prove  that  for  any  function  tp  periodic  on  T>t  the  following  weak  form  is 
satisfied 


Q  (/,  /)  <p{v)  dv  =  - 


/  /  B(x,y)8(x-y) 

I  V>t  ^  J  vGT>t  J  xGBR  J y£BR 

f(v  +  x  +  y)f(v)  [p{v  +  y)  +  (p(v  +  x)  —  <p(v  +  x  +  y)  —  <p{v)\  dv  dx  dy.  (44) 


About  the  conservation  properties  one  can  shows  that  if  /  has  compact  support  in¬ 
cluded  in  Br  with  T  >  (3  + y/2)/?/2  (no  aliasing  condition,  see  Pareschi  and  Russo  (2000b) 
for  a  detailed  discussion),  then  no  unphysical  collisions  occur  and  thus  mass,  momentum 
and  energy  are  preserved.  Obviously  this  compactness  is  not  preserved  with  time  since 
the  collision  operator  spreads  the  support  of  /  by  a  factor  y/2.  In  the  rest  of  the  paper 
we  will  focus  on  the  periodized  truncation  QR. 
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Figure  5:  Periodization  in  2D 

2.2  Spectral  methods 

Now  we  use  the  representation  QR  to  derive  new  spectral  methods.  The  spectral  methods 
for  kinetic  equations  originated  in  the  works  of  Pareschi  and  Perthame  (1996),  Pareschi 
and  Russo  (2000b),  and  were  further  developed  in  Pareschi  and  Russo  (2000c)  and  Filbet 
and  Russo  (2003).  Before  they  had  a  long  history  in  fluid  mechanics,  see  Canuto  et  al. 
(1988). 

To  simplify  notations  let  us  take  T  =  n.  Hereafter  we  use  just  one  index  to  denote 
the  d-dimensional  sums  of  integers. 

The  approximate  function  /jy  is  represented  as  the  truncated  Fourier  series 

N 

!n(v)  =  Y, 
k=—N 

The  spectral  equation  is  the  projection  of  the  collision  equation  in  PjV,  the  (2 N  +  l)d- 
dimensional  vector  space  of  trigonometric  polynomials  of  degree  at  most  N  in  each  direc¬ 
tion,  i.e 

=  'P nQRUn ,  In) 

where  Vn  denotes  the  orthogonal  projection  on  PA  in  L2(T>n).  A  straightforward  compu¬ 
tation  leads  to  the  following  set  of  ordinary  differential  equations  on  the  Fourier  coefficients 

N 

fk(t)=  ^2  k  =  -N,...,N  (45) 

l,m=—N 

l-\-m=k 

where  (3(1,  m)  are  the  so-called  kernel  modes,  given  by 

P(l,  m)  —  f  I  B(x,  y )  5(x  ■  y )  [eil'x  eim  y  -  einHx+y)]  dx  dy. 

J x£Br  J y&Bn 
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The  kernel  modes  can  be  written  as 

(3(1,  m)  =  (3(1 ,  m)  —  (3(m,  m) 


where 

(3(1,  m)  —  f  I  B(x,y)8(x-y)eil  xeim  ydxdy. 

Jx£Br  J V&Br 

Therefore  in  the  sequel  we  shall  focus  on  (3,  and  one  easily  checks  that  (3(1,  m)  depends 
only  on  |/|,  \m\  and  1 1  ■  m\. 

Finally  let  us  compare  the  new  kernel  modes  with  the  ones  in  Pareschi  and  Russo 
(2000b).  The  usual  kernel  modes  written  in  the  x  and  y  variables  reads 

Vualfi.  m)  =  f  f  B(x,  y)  S(x  ■  y)  x{k+„(gs(  [e“-  e‘”*  »  -  dx  dy. 

JxGBr  j V&Br 

Thus  the  usual  representation  contains  a  strong  coupling  between  x  and  y  which  makes 
it  very  hard  the  construction  of  fast  algorithms. 

2.3  Discrete-velocity  models 

The  representation  QR  of  this  section  can  also  be  used  to  derive  fast  solvers  for  discrete 
velocity  models  (DVM).  ffistorically  these  methods  were  among  the  first  deterministic 
methods  for  discretizing  the  Boltzmann  equation  in  velocity  space.  The  discretization 
is  built  starting  from  physical  rather  then  numerical  considerations.  We  assume  the  gas 
particles  can  attain  only  a  finite  set  of  velocities 

VN  =  {vi,  v2,  v3, ...,  vN},  Vi  G  JR3. 

Any  DVM  can  be  written  as  a  product  quadrature  formula  for  (3)  in  the  general  form 

d,~  V  rfj  [ftf,  -  f,fj] , 

j,k,l  G 

where  Di  denotes  the  discrete  Boltzmann  collision  operator  and  the  integer  indexes  refer 
to  the  points  in  the  computational  grid.  In  order  to  keep  conservations  the  coefficients 
rfj  are  defined  by 

rfj  =  1(*  +  j  -  k - /)  1(|*|2  +  \j\2  -  \k\2  -  |/|2) B(\k j\) 

where  1  denotes  the  function  on  Z  defined  by  l(^)  =  1  if  z  —  0  and  0  elsewhere,  and 
>  0  are  the  weights  of  the  quadrature  formula,  which  characterize  the  different  DVM. 
B  >  0  is  the  discrete  collision  kernel.  One  can  check  on  this  formulation  that  the  scheme 
satisfies  the  usual  conservation  laws  and  entropy  inequality  (see  Platkowski  and  Illner 
(1988)  and  the  references  therein). 

We  can  write  at  the  discrete  level  the  same  representation  as  in  the  continuous  case 

Di  ^  ^  f  k,l  \fi+kfi+l  fifi+k+l\ 

k,iezd 
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Figure  6:  Sketch  of  a  planar  model  based  on  a  cartesian  grid.  Note  that  in  general  few 
grid  points  will  belong  to  the  collision  circle. 


with 


Tfc,;  —  1  (k 


()1F0T  w' 


This  is  coherent  with  the  DVM  obtained  by  quadrature  starting  from  the  Carleman  rep¬ 
resentation  in  Panferov  and  Heintz  (2002). 

Now  again  when  one  is  interested  to  compute  the  DVM  in  a  bounded  domain  there 

-f-y* 

are  two  possibilities.  First  as  in  the  case  of  QhL  one  can  force  the  discrete  velocities  to 
stay  in  a  box,  which  yields  for  i  =  —  N, . . . ,  N  (again  using  the  one  index  notation  for 
d-dimensional  sums) 


^  ^  f  k,l  \fi+kfi+l  fifi+k+l\  • 


k,l 

—  N<  i+fc,  i-\-l ,  i-\-k-\-l<N 


This  new  discrete  operator  is  completely  conservative  but  the  collision  kernel  is  not  in¬ 
variant  anymore  according  to  i,  which  breaks  the  convolution  properties. 

The  other  possibility  is  to  periodize  the  function  /  over  the  box  and  truncate  the  sum 
in  k  and  l.  It  yields  for  a  given  truncation  parameter  N  6  N 

D?  =  Y.  ^[fi+kfi+i-fifi+k+i],  (46) 

-N<k,l<N 


for  any  i  =  —  N . . .  N . 

It  is  easy  to  see  that  DN  satisfies  exactly  a  discrete  weak  form  and  conservation  proper¬ 
ties  similar  to  QR.  Moreover  one  can  derive  the  following  consistency  result  from  (Panferov 
and  Heintz,  2002,  Theorem  3)  in  the  case  of  hard  spheres  collision  kernel 
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2.4  Fast  spectral  algorithms 


Theorem  1  Assume  that  f,g  e  Cfc(K3)  (k  >  1)  with  compact  support  Bs ■  The  uniform 
grid  of  step  h  is  constructed  on  the  boxVj-  with  the  no  aliasing  condition  T  >  (3+v/2)S'/2. 
Then  for  N  =  [S/h]  (where  [•]  denotes  the  integer  value)  and  h  >  0  sufficiently  small, 


\\Q(s,  f)  ~  (g,  f)\\L°°(zh)  <  Chr 

where  D(f  is  the  DVM  operator  defined  in  (46)  (for  the  precise  quadrature  weights  derived 
in  Panferov  and  Heintz  (2002))  on  the  grid  above-mentioned,  and  =  f(ih).  Here 
r  =  k/(k  +  3)  and  the  constant  C  is  independent  on  h. 


Remark  2  As  can  be  seen  from  Theorem  1,  the  periodized  DVM  presented  in  this  subsec¬ 
tion  is  expected  to  have  a  quite  low  accuracy.  On  the  contrary  the  spectral  method  will  be 
proven  to  be  spectrally  accurate,  i.e.  of  infinite  order  for  smooth  solutions.  Nevertheless 
this  periodized  DVM  has  some  interesting  features  compared  to  the  spectral  method.  In¬ 
deed,  one  can  prove  that  if  the  quadrature  weights  Wf-j  are  non-negative,  then  the  scheme 
is  stable  in  the  sense  that  if  one  starts  from  a  non-negative  initial  data,  then  the  solution 
remains  non-negative  and  has  thus  constant  L 1  norm.  Concerning  the  spectral  method 
we  refer  to  Filbet  and  Mouhot  (2011)  for  an  analysis  of  its  stability  and  convergence 
properties. 


2.4  Fast  spectral  algorithms 

As  soon  as  one  is  searching  for  fast  deterministic  algorithms  for  the  collision  operator,  i.e 
algorithm  with  a  cost  lower  than  0(N2d+£ )  (which  is  the  cost  of  a  usual  discrete  velocity 
model,  with  typically  £  =  1),  one  has  to  find  some  way  to  compute  the  collision  operator 
without  going  through  all  the  couples  of  collision  points  during  the  computation.  This 
leads  naturally  to  search  for  some  convolution  structure  (discrete  or  continuous)  in  the 
operator.  Unfortunately,  as  discussed  in  the  previous  sections,  this  is  rather  contradictory 
with  the  search  for  a  conservative  scheme  in  a  bounded  domain,  since  the  boundary 
condition  needed  to  prevent  for  the  outgoing  or  ingoing  collisions  breaks  the  invariance. 

Here  we  search  for  a  convolution  structure  in  the  equations  (45).  The  aim  is  to 
approximate  each  (3(1,  m)  by  a  sum 


A 

(3(1, m )  ~  ^ otp(l)otv(m). 

p= i 


This  gives  a  sum  of  A  discrete  convolutions  and  so  the  algorithm  can  be  computed  in 
0(A  Nd  log2  N)  operations  by  means  of  standard  FFT  techniques  (Canuto  et  ah,  1988; 
Cooley  and  Tukey,  1965).  Obviously  this  is  equivalent  to  obtain  such  a  decomposition 
on  (3.  To  this  purpose  we  shall  use  a  further  approximated  collision  operator  where  the 
number  of  possible  directions  of  collision  is  reduced  to  a  finite  set. 
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2.4.1  A  semi-discrete  collision  operator 

We  write  x  and  y  in  spherical  coordinates 

QR(fJ)(v)=  J  [  [  8{e-e')dede' 

4  deG Sd_1  Je'eSd_1 

f  [  Pd~ 2  (p')d“2  B(p,  p')  [f(v  +  p'e')f{y  +  pe)  -  /(u  +  pe  +  p'e')/(v)]  dp  dp' 

J-R  J-R 

(47) 

Let  us  take  A  a  discrete  set  of  orthogonal  couples  of  unit  vectors  (e,e7),  which  is  even: 
(e,  e')  G  A  implies  that  (— e,  e'),  (e,  —  e')  and  (— e,  — e')  belong  to  A  (this  property  on  the 
set  A  is  required  to  preserve  the  conservation  properties  of  the  operator).  Now  we  define 
Qx  to  be 


QRAU,f)(v)  =  ^ 


pR  pR 


'(e,e')e*4  I  J-R  J-R 


pd-2  {p,)d-2  [/(u  +  p'e7)/^  +  pe)  — 


/(u  +  pe  +  p'e')/(u)]  dpdp7|  dA 


where  dA  denotes  a  discrete  measure  on  A  which  is  also  even  in  the  sense  that  dA(e,  e')  = 
dA(— e,e')  =  dA(e,  — e')  =  dA(— e,  —  e').  Using  again  translation  change  of  variable  on  v 
by  pe,  pV  and  pe  +  p'e'  and  the  symmetries  of  the  set  A  one  can  easily  derive  the  following 
weak  form  on  Q^.  For  any  function  <p  periodic  on  T)?, 


r*R  pR 


fVT 


QHrA(fJ)<p(v)dv  =  - 

lo 


>veVT  J(e,e')&A  J-R  J-R 


Pd~2  (p'j  -2  B(p,  p') 


f(v  +  pe  +  p'e')f{v )  [<p{v  +  p' e')  +  <p(v  +  pe)  —  <p(u  +  pe  +  p'e')  —  <p(u)]  dp  dp'  dAdv. 


This  immediately  gives  the  same  conservations  properties  as  Qr. 


2.4.2  Expansion  of  the  kernel  modes 

We  make  the  decoupling  assumption  that 

B(x,y)  =  a(\x\)b(\y\).  (48) 

This  assumption  is  obviously  satisfied  if  B  is  constant.  This  is  the  case  of  Maxwellian 
molecules  in  dimension  two,  and  hard  spheres  in  dimension  three  (the  most  relevant  kernel 
for  applications).  Extensions  to  more  general  interactions  are  discussed  in  Mouhot  and 
Pareschi  (2006). 

First  let  us  deal  with  dimension  2  with  B  —  1  to  explain  the  method.  Here  we  write 
x  and  y  in  spherical  coordinates  x  =  pe  and  y  =  pV  to  get 


=  jf  [  d(e-e') 

r  rR  i 

/  eip{l'e)  dp 

''CD 

A 

de  de'. 

4  dees1  de'GS1 

J-R 

L  J-R  J 
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Let  us  denote  by 


rR 


JR 


S  = 


eips  dp , 


l-R 


for  sGl.  It  is  easy  to  see  that  <j)2R  is  even  and  we  can  give  the  explicit  formula 


JR 


s)  —  2  R  Sinc(-Rs) 


with  Sinc(d)  =  (sin0)/0. 
Thus  we  have 


'ees1  Je'eS1 


S(e  •  e7)  (j)R{l  •  e)  4>R(m  •  e ')  dede' 


and  thanks  to  the  parity  property  of  <fi2R  we  can  adopt  the  following  periodic  parametriza- 
tion 

/»7T 

P(l,  m)  =  /  <j>\{l  ■  eg)  4>2R{m  ■  ee+n/2)  dO. 

J  o 

The  function  6  — »  <f>2R(l  ■  eg)  4>2R[m  ■  eg+n/2)  is  periodic  on  [0, 7r]  and  thus  the  rectangular 
quadrature  rule  is  of  infinite  order  and  optimal.  A  regular  discretization  of  M  equally 
spaced  points  thus  gives 

M—l 

P(l,m)  =  jj  «p(0«pM 

p= o 

with 

«p(0  =  <Pr(1  ■  e9P)r  <x'p(m)  =  (f)R(m  ■  egp+n/2) 

where  9p  =  np/M. 

More  generally  under  the  decoupling  assumption  (48)  on  B ,  we  get  the  following 
decomposition  formula 

M—l 

P(l,m)  =  J] 

p= o 

ap(m)  —  (t>2R,b(m  '  e0P+i r/2) 

<i>R,b(s)  =  [  b{p‘)  etf/a  dp’ 

J-R 

with  9p  =  tt p/M. 


where 

and 


(Xp{l)  0_R,a(^  ‘  e9P)i 


r-  R 


JR,ci 


S  = 


a(p )  elps  dp, 


Remark  3  In  the  symmetric  case  a  =  b  (for  instance  for  hard  spheres)  it  is  possible  to 
parametrize  (3(1,  m)  as 

fir/2 

P(l,  m)  =  2  /  (fR}a(l  ■  eg)  4>\a(m  ■  e0+7r/2)  dO 
Jo 

and  the  function  9  — »  (f2R  a(l  ■  eg)  <fRa(m  ■  ep+7r/2)  is  periodic  on  [0, 7r/2] .  Thus  the  decompo¬ 
sition  can  be  obtained  by  applying  the  rectangular  rule  on  this  interval.  At  the  numerical 
level  it  yields  a  reduction  of  the  cost  by  a  factor  2. 
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Now  let  us  deal  with  dimension  d  —  3  with  B  satisfying  the  decoupling  assump¬ 
tion  (48).  First  we  change  to  the  spherical  coordinates 


P(!,m)  =  - 


Ue S2  Je'es2 


5(e  ■  e') 


rR 


pa(p)  e dp 


l-R 


rR 


p'  b(p')  dp' 


I  —R 


de  de' 


and  then  we  integrate  first  e'  on  the  intersection  of  the  unit  sphere  with  the  plane  eL , 


=  j  I  •  e) 

4  lee §2 


L  J  e/GS2ne-L 


<f)3R  h(rn  ■  e')  de' 


de 


where 


rR 


JR,a 


S  = 


pa(p)  eips  dp. 


I-R 


Thus  we  get  the  following  decoupling  formula  with  two  degrees  of  freedom 


P(l,m)=  /  <j>^a(l-e)^b(lle±(m))de 

Jee  Si 


where  denotes  the  half-sphere  and 

ipR,b( ne±(m))  =  /  sin6*0Rift(|ne±(m)|  cos6>)  dO, 

Jo 

(this  formula  can  be  derived  performing  the  change  of  variable  de'  =  sin  9  dQ  dip  with  the 
basis  (e,  u  =  ne±(m)/|IIe±(m)|,  e  x  u)). 

Again  in  the  particular  case  where  B  —  1  (hard  spheres  model),  we  can  compute 
explicitly  the  functions  <j)R  (in  this  case  a  =  b  =  1), 

</>|(s)  =  R2  [2Sinc (Rs)  -  Sinc2(i?s/2)]  ,  i/j3r(s)  =  2  R2  Sinc2(Rs/2). 

Now  the  function  e  — >■  (ftR  a(l  ■  e)  i/)R  b (lle±  (m))  is  periodic  on  §+  and  so  the  rectangular 
rule  is  of  infinite  order  and  optimal.  Taking  a  spherical  parametrization  ( 9 ,  p)  of  e  G  S>^ 
and  uniform  grids  of  respective  size  Mi  and  M2  for  9  and  p  we  get 


2  Mi,M2 


J(l-in)  1M/  oip,q(l)a'pq(m) 


p,q= o 


where 

and 

and 


0  =  ■  e(eP)V,)),  dp>q{m)  =  ip%b(n e±^(m)) 


rR 


<t>R,a(s)  =  /  pa(p)elpsdp,  ipR,b(s)=  sin 9  (j)\b(s  cos 9)  d9 


l-R 


fa  i  _  fPn  Q^\ 

(6ln  Vq)  \Mp  mJ' 


From  now  on  we  shall  consider  this  expansion  with  M  =  M)  =  M2  to  avoid  anisotropy  in 
the  computational  grid. 
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Remark  4 

For  any  dimension,  we  can  construct  as  above  an  approximated  collision  operator 
qR,Am  Wllh 

Am  =  {(e,  e')  e  x  |  e  G  $££,  eAe1!!  S^1} 

where  denotes  a  uniform  angular  discretization  of  the  half  sphere  with  M  points  in 
each  angular  coordinate  (the  other  half  sphere  is  obtained  by  parity).  Let  us  remark  that 
this  discretization  contains  exactly  M^1  points.  From  now  on  we  shall  denote 


QRM  =  Q 


_Md_1 

RAm  =  ^  q 

P=  1 


RM 

P 


2.4.3  Spectral  accuracy 

In  this  paragraph  we  are  interested  in  computing  the  accuracy  of  the  scheme  according 
to  the  three  parameters  N  (the  number  of  modes),  R  (the  truncation  parameter),  and 
M  (the  number  of  angular  directions  for  each  angular  coordinate).  Instead  of  looking  at 
the  error  on  each  kernel  mode  it  is  more  convenient  to  look  at  the  error  on  the  global 
operator.  Here  the  Lebesgue  spaces  Lp,  p  =  1 . . .  +  oo,  and  the  periodic  Sobolev  spaces 
H!f,  k  —  0  . . .  +  oo  refer  to  Vw. 

So  in  order  to  give  a  consistency  result,  the  first  step  will  be  to  prove  a  consistency 
result  for  the  approximation  of  QR  by  QR,M  (see  Mouhot  and  Pareschi  (2006)  for  details). 


Lemma  1  The  error  on  the  approximation  of  the  collision  operator  is  spectrally  small, 
i.e  for  all  k  >  d  —  1  such  that  f  e 


\\Qa(s,f)  -  <2r'm(j./)IU»  <  Ci 


Mk 


For  the  second  step  we  shall  use  the  consistency  result  (Pareschi  and  Russo,  2000b, 
Corollary  5.4)  on  the  operator  QR,  which  we  quote  here  for  the  sake  of  clarity. 

Lemma  2  For  all  k  e  N  such  that  f  e 


\\Qr(U)-VnQrUn.!A\W  <  ff  (ll/llwf  +  WQVnJAWhi 

Combining  these  two  results,  one  gets  the  following  consistency  result 
Theorem  2  For  all  k  >  d  —  1  such  that  f  e  Ft*  (D^) 

Rk\\fN\\lk  r  / 

\\QR(fJ)  -VNQR’M(fN,fN)\\L2  <  Cy  Mk  p  (h/IIhj  +  \\QR(fNjN)\\H! 

Now  let  us  focus  briefly  on  the  macroscopic  quantities.  First  with  Lemma  1  at  hand 
one  can  establish  the  estimate 


\\QR-M(9J)\\li<C 


Hd 


8-28 


RTO-EN-AVT-194 


2.5  Fast  DVM’s  algorithms 


2  FAST  BOLTZMANN  SOLVERS 


for  a  constant  uniform  in  M.  Then  following  the  method  of  (Pareschi  and  Russo,  2000b, 
Remark  5.4)  and  using  this  estimate  we  obtain  the  following  spectral  accuracy  result 

-  (PNQR'M(fN,fN),v)\L1  <  IMU=  (ll/ll„;«  +  II«s-mUWaOIIh}) 

where  tp  can  be  replaced  by  v,  |u|2.  Indeed  there  is  no  need  to  compare  the  momenta 
of  VnQr,m (/tv,  In)  with  those  of  QR(f,f )  since  QR,M  is  also  conservative,  and  so  they 
can  be  compared  directly  to  those  of  QR,M .  Thus  the  error  on  momentum  and  energy  is 
independent  on  M  and  is  spectrally  small  according  to  N  even  for  very  small  value  of  the 
parameter  M . 

2.4.4  Implementation  aspects 

The  final  spectral  scheme  depends  on  the  three  parameters  N,  R ,  and  M .  The  only 
conditions  on  these  parameters  is  the  no-aliasing  condition  that  relates  R  and  the  size  of 
the  box  T  (here  n).  A  detailed  study  of  the  influence  of  the  choices  of  N  and  R  has  been 
done  in  Pareschi  and  Russo  (2000b).  Here  we  are  interested  only  in  the  influence  of  M 
over  the  computations,  since  M  controls  the  computations  speed-up. 

The  method  of  the  previous  subsections  yields  a  decomposition  of  the  collision  opera¬ 
tor,  which  after  projection  on  PjV  gives  the  following  decomposition 

Md_1 

VnQr’m  =  Y,  VnQr'm.  (49) 

p= i 

Each  VnQr’m  can  be  computed  with  a  cost  0(Nd  log2  N).  Thus  for  a  general  choice  of 
M  and  N  we  obtain  the  cost  0(A/Id~1Nd  log2  N).  The  decomposition  (49)  is  completely 
parallelizable  and  thus  the  cost  can  be  strongly  reduced  on  a  parallel  machine  (formally 
up  to  0(Nd  log2  N)).  One  just  has  to  make  independent  computations  for  the  Md  l  terms 
of  the  decomposition. 

Moreover  the  formula  of  decomposition  is  naturally  adaptive  (that  is  the  number  M 
can  be  made  space  dependent),  which  can  be  quite  useful  in  the  inhomogeneous  setting, 
where  some  regions  deserve  less  accuracy  than  others.  Since  it  relies  on  the  rectangular 
formula,  whose  adaptivity  property  is  well  known,  one  can  easily  double  the  number  of 
directions  M  if  needed,  without  computing  again  those  points  already  computed. 

Finally  the  decomposition  can  be  also  interesting  from  the  storage  viewpoint,  as  the 
classical  spectral  method  requires  the  storage  of  a  Nd  x  Nd  matrix  whereas  the  fast 
method  requires  the  storage  of  2Md~1  vectors  of  size  Nd.  In  dimension  2  the  classical 
method  requires  a  storage  of  order  0(N 4)  and  the  fast  method  requires  a  storage  of  order 
0(MN2).  In  dimension  3  the  classical  method  requires  a  storage  of  order  0(N 4)  (thanks 
to  the  symmetries  of  the  matrix  of  kernel  modes,  see  Pareschi  and  Russo  (2000b)),  and 
the  fast  method  requires  a  storage  of  order  0(M2  N3). 

2.5  Fast  DVM’s  algorithms 

The  fast  algorithms  developed  for  the  spectral  method  can  be  in  fact  extended  to  the 
periodized  DVM  method.  The  method  that  originates  is  in  some  sense  related  to  the 
direct  FFT  approach  proposed  in  Bobylev  and  Rjasanow  (1997,  2000,  1999). 
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2.5.1  Principle  of  the  method:  a  pseudo-spectral  viewpoint 

We  start  from  the  periodized  DVM  in  [|  —  N,  N\]d  with  representation  (46)  and  as  in  the 
continuous  case  we  set,  for  —V  <  k,l  <  N, 


B(\k\,\l\) 


mm 

|  k  +  l\d~ 2 


■)d—l 


B 


1*1 


vfF+lf 


,VW+ 


+ 


With  this  notation 

f k,i  =  Hk  ■  l)B(\k\,  | l\)wKh 

and  thus  the  DVM  becomes 


fi=  E  1(fc  '  0  B(\k\,\l\)  wKl  [ fi+kfi+i  -  fifi+k+i]  ■ 

- N<k,l<N 


Now  we  transform  this  set  of  ordinary  differential  equations  into  a  new  one  using  the 
involution  transformation  of  the  discrete  Fourier  transform  on  the  vector  ( fi)-N<i<N ■ 
This  involution  reads 

2N  N 

h  =  .,v  !  E^e-^)’  fi  =  E  ~hei (*) 

1=0  I=—N 


Zl 7T  -K 

where  exyk )  denotes  e  2N+ 1  ,  and  thus  the  set  of  differential  equations  becomes 


N 


2  N 


f,  =  £ 

K,L=—N 


2N  +  1 


E  eR+L 


-IW 


i= 0 


E  l(fc-i)  B(\k\,  \l\)wKi  ( eK(k)eL(l)  -  eL(k  +  l ) 

-N<k,l<N 


Ik  /l 


for  —  N  <  I  <  N.  We  have  the  following  identity 

1  2  N 

— —  E  ^K+L-i(i)  =  1(AT  +  L  -  /) 

^  i=0 

and  so  the  set  of  equations  is 


N 

?,=  £  P(K,L)fKfL 

K+L=I 

K,L=-N 


with 

m,  L)  =  E  ■  0  Kl)  wk,i  [eK(k)eL(l)  -  eL(k  +  /)]  =  /3(W,  L)  -  /3(L,  L) 

-N<k,l<N 

where 

/3(A',  L)  —  E  !(*  ■  0  I)  eK(k)eL(l)- 

- N<k,l<N 


8-30 


RTO-EN-AVT-194 


2.5  Fast  DVM’s  algorithms 


2  FAST  BOLTZMANN  SOLVERS 


Let  us  first  remark  that  this  new  formulation  allows  to  reduce  the  usual  cost  of  compu¬ 
tation  of  a  DVM  exactly  to  0(N2d)  (as  with  the  usual  spectral  method).  Note  however 
that  the  (2 A  +  l)d  x  (2 A  +  l)d  matrix  of  coefficients  (/3(K,  L))k,l  has  to  be  computed 
and  stored  first,  thus  the  storage  requirements  are  larger  with  respect  to  usual  DVM. 
Now  the  aim  is  to  give  an  expansion  of  /3(K,  L)  of  the  form 

M 

Pk,l  ^J^ap(K)  a'p(L) 
p= i 

to  get  a  lower  cost  by  the  use  of  discrete  convolution. 


2.5.2  Expansion  of  the  discrete  kernel  modes 

We  make  a  decoupling  assumption  as  in  the  spectral  case 

B(\k\,\l\)wkil  =  a(k)b(l). 

Remark  5  Note  that  the  DVM  constructed  by  quadrature,  in  dimension  3  for  hard  spheres, 
in  Panferov  and  Heintz  (2002)  satisfies  this  decoupling  assumption  with  a(k )  =  hA / gcd(ki,  k2 ,  k3) 
and  b(l)  =  1  (see  (Panferov  mid  Heintz,  2002,  Formula  (2.8))),  and  gcd(ki,  k2,  k%)  denotes 
the  greater  common  divisor  of  the  three  integers. 


The  difference  here  with  the  spectral  method,  which  is  a  continuous  numerical  method, 
is  that  we  have  to  enumerate  the  set  of  {— N  <  k,l  <  N  \  k  HI}.  This  motivates  for  a 
detailed  study  of  the  number  of  lines  passing  through  0  and  another  point  in  the  grid.  To 
this  purpose  let  us  introduce  the  Farey  series  and  a  new  parameter  0  <  A  <  N  for  the 
size  of  the  grid  used  to  compute  the  number  of  directions.  The  usual  Farey  serie  is 


Fh  =  {(p>9)  e  [\°,N\?  I  0  <  p  <  q  <  N  and  gcd (p,q)  =  1} 

where  gcd(p,  q)  denotes  again  the  greater  common  divisor  of  the  two  integers  (more  details 
can  be  found  in  Hardy  and  Wright  (1979)).  It  is  straightforward  to  see  that  the  number 
of  lines  Ah  passing  through  0  in  the  grid  [|  —  A,  A|]2  is  Ak  =  4  |V(  |.  In  fact  symmetries 
often  allow  to  reduce  the  number  of  directions  needed.  Similarly  one  can  define 


=  { (p,  q,  r)  G  [|  0,  N\Y’  |  0  <p<q<r<N  and  gcd(p,  q,  r)  —  l} 

and  the  number  of  lines  A k  passing  through  0  in  the  grid  [|  —  iV,iV|]3  is  A k  =  16 
The  exponents  of  the  Farey  series  refer  to  the  dimension  of  the  space  of  lines  (which  is 
d  —  1).  Now  let  us  estimate  the  cardinal  of  J-T  and  (see  Mouhot  and  Pareschi  (2011) 
for  details). 


Lemma  3  The  Farey  series  in  dimension  d  —  2  and  d  =  3  satisfy  the  following  asymptotic 
behavior 


J7- 

'rN 

T2- 

■J  N 


N2 

2C(2) 

N3 

403) 


+  0(N  log  A”) 
+  0(N2) 


3  N2 


7T- 


+  0(N  log  A) 


where  (  denotes  the  usual  zeta  function. 
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2.6  Numerical  results 


Now  one  can  deduce  the  following  decomposition  of  the  kernel  modes 
D(K,L)  =  J2  ViWWWI'DMUeaO 

-N<k,l<N 


-  V 

Y  a(\k\)  eK(k) 

E  KI'IK  (0 

e€AN 

-  k£e7j,  —N<k<N 

-lee-L^NKlKN 

with  equality  if  TV  =  TV.  Here  Afj  denotes  the  set  of  primal  representants  of  directions  of 
lines  in  [|  —  TV,  TV|]  passing  through  0.  After  indexing  this  set,  which  has  cardinal  A^,  one 
gets 

K 

(5°) 

p= i 

with 

<*p(k)=  Y  a(\k\)eK{k),  oi'p(L)  =  Y  b(\l\ )ei(0- 

keepZ,-N<k<N  lee£,-N<l<N 


2.5.3  Implementation  aspects 


The  method  yields  a  decomposition  of  the  discrete  collision  operator 


A% 


Dn  ~  Dn’R  =  E  j~^N,N,p 
P=  1 


with  equality  if  TV  —  TV.  Each  DN,N,p(f,  f)  is  defined  by  the  p-th  term  of  the  decomposition 
of  the  kernel  modes  (50).  Each  term  DN,N,P  of  the  sum  is  a  discrete  convolution  operator 
when  it  is  written  in  Fourier  space. 

Thus  one  can  see  that  even  if  we  take  TV  =  TV  =  N,  i.e  we  take  all  possible  directions 
in  the  grid  [|  —  N,  iV|]d,  we  get  the  computational  cost  0(N2d  log2  N)  which  is  better  than 
the  usual  cost  of  the  DVM,  0(N2d+1 )  (but  slightly  worse  than  the  cost  0(N2d)  obtained 
by  solving  directly  the  pseudo-spectral  scheme,  thanks  to  a  bigger  storage  requirement). 

More  generally  for  a  choice  of  TV  <  N  we  obtain  the  cost  0(NdNd  log2  TV).  The  same 
remarks  we  did  for  the  fast  spectral  algorithms  about  the  parallelization  and  adaptivity 
(and  storage  interest)  of  the  method  hold  true  in  this  case:  a  parallel  algorithm  could 
reduce  the  computational  cost  up  to  0(Nd  log2  TV). 

Moreover  we  expect  that  for  DVM  one  can  strongly  reduce  the  parameter  TV  in  order 
to  improve  the  cost  of  the  scheme,  without  damaging  the  accuracy  of  the  scheme.  The 
justification  for  this  is  the  low  accuracy  of  the  method  (the  reduction  of  the  number  of 
direction  has  a  small  effect  on  the  already  poor  accuracy  of  the  scheme). 


2.6  Numerical  results 

In  this  section  we  will  present  several  numerical  results  for  the  space  homogeneous  equation 
which  show  the  improvement  of  the  fast  spectral  algorithms  with  respect  to  the  classical 
spectral  methods.  The  time  discretization  is  performed  by  standard  explicit  Runge-Kutta 
methods. 
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2.6.1  2D  Maxwell  molecules 
Comparison  to  exact  solutions 

We  consider  2D  pseudo-Maxwell  molecules  (i.e.,  the  VHS  model  with  7  =  0).  In  this  case 
we  have  an  exact  solution  given  by 


f{t,v)  = 


exp(— u2/2  5) 
2vr  52 


25-1  + 


1-5 

25 


with  5  =  1— exp(— 1/8)/2,  which  corresponds  to  the  well  known  “BKW”  solution  (Bobylev, 
1975).  This  test  is  performed  to  check  spectral  accuracy,  by  comparing  the  error  at  a 
given  time,  when  using  nv  =  8,  16  and  32  Fourier  modes  for  each  coordinate.  We  present 
the  results  obtained  by  the  classical  spectral  method  and  the  fast  spectral  method  with 
different  numbers  of  discrete  angles. 

Figure  7  shows  the  relative  L°°,  L  x ,  and  L2  norms  of  the  difference  between  the  nu¬ 
merical  and  the  exact  solution,  as  a  function  of  time.  We  refer  to  Filbet  et  al.  (2006)  for 
a  more  detailed  discussion  about  the  different  source  of  error. 

Concerning  the  comparison  between  the  classical  and  fast  spectral  methods,  we  observe 
that  for  a  fixed  value  of  nv ,  the  numerical  error  of  the  classical  spectral  method  and  of 
the  fast  algorithm  is  of  the  same  order.  Moreover,  the  influence  of  the  number  of  discrete 
angles  is  very  weak.  In  Table  1,  we  give  a  quantitative  comparison  of  the  numerical  error 
Si  at  time  Tend  =  1.  We  can  also  observe  the  spectral  accuracy  for  the  classical  and 
fast  methods:  the  order  of  accuracy  is  about  3  between  8  and  16  grid  points,  whereas  it 
becomes  7  between  16  to  32  points. 


Efficiency  and  accuracy 


Now,  we  still  consider  2D  pseudo-Maxwell  molecules  (i.e.,  7  =  0) 
datum 


/( 0,u) 


+  exp 


^  K>  +  t>o|L> 


with  the  following  initial 


v  e  M2, 


where  vq  =  (1,2).  In  this  case,  we  do  not  know  the  exact  solution  but  we  want  to  study 
the  influence  of  the  number  of  discrete  angles  on  a  non-isotropic  solution.  Thus,  this  test 
is  used  to  check  the  energy  conservation  and  the  evolution  of  high-order  moments  of  the 
solution.  Simulations  are  performed  with  nv= 16,  32  and  64  points. 


Number  of 
points 

Classical 

spectral 

Fast  spectral 
with  M  =  4 

Fast  spectral 
with  M  —  6 

Fast  spectral 
with  M  —  8 

8 

0.02013 

0.02778 

0.02129 

0.02112 

16 

0.00204 

0.00329 

0.00238 

0.00224 

32 

1.405E-5 

2.228E-5 

1.861E-5 

1.772E-5 

Table  1:  Comparison  of  the  L 1  error  in  2 D  between  the  classical  spectral  method  and  the 
fast  spectral  method  with  different  numbers  of  discrete  angles  and  with  a  second-order 
Runge-Kutta  time  discretization  at  time  Tend  =  1. 
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2  4  6  8  10  12  14  2  4  6  8  10  12  14 


Figure  7:  Evolution  of  the  numerical  L 1  and  L°°  relative  error  of 

In  Figure  8  the  relaxation  of  the  entropy  and  the  temperature  components  for  the  fast 
and  classical  spectral  methods  is  shown.  Finally,  we  plot  in  Figure  9  the  time  evolution 
of  high-order  moments  of  given  in  discrete  form  by 

N 

Mk(t)  =  Av2  Y  \vi\k  fN(t,vi). 

l=—N 

High-order  moments  give  information  on  the  accuracy  of  the  approximate  distribution 
function  tail.  Once  again,  we  observe  that  the  number  of  angles  does  not  affect  the 
results  even  if  the  solution  is  non-isotropic. 


Figure  8:  Relaxation  of  the  entropy  and  the  temperature  components  for  the  fast  and 
classical  spectral  methods  with  respect  to  the  number  of  modes  per  direction  nv  and  the 
length  box  T. 
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Figure  9:  Time  evolution  of  the  variations  of  high  order  normalized  moments  M. 4, 
and  of  f(t,v )  for  the  fast  and  classical  spectral  methods  with  respect  to  the  number 
of  modes  per  direction  nv  and  the  length  box  T. 


2.6.2  3D  Hard  Spheres 


In  this  section  we  consider  the  3D  Hard  Sphere  molecules  (HS)  model.  The  initial  condi¬ 
tion  is  chosen  as  the  sum  of  two  Gaussians 


f(v,  0) 


1 

2(2vru2) 


exp 


+  exp 


with  (j  =  l  and  Vq  =  (2, 1,  0).  The  final  time  of  the  simulation  is  Tenc[  =  3  and  corresponds 
approximatively  to  the  time  for  which  the  steady  state  of  the  solution  is  reached.  The 
time  step  is  At  =  0.1  and  the  length  box  is  taken  as  T  =  12  when  nv  =  16  and  T  =  15 
when  nv  =  32. 

This  test  is  used  to  check  the  evolution  of  moments  and  particularly  the  stress  tensor 
Pij,  i,j  =  1,  •  •  •  ,3  defined  as 


P 

1 


f{v){vi  -  Ui )  (Vj  -  Uj )  dv,  ( i,j )  e  {1,  2,  3}2, 


where  (iq)j  are  the  components  of  the  mean  velocity.  In  Figure  10,  we  propose  the  evolu¬ 
tion  of  the  temperature  for  the  two  methods  using  32  grid  points  in  each  direction.  The 
solution  is  also  compared  with  the  solution  obtained  from  a  standard  Direct  Simulation 
Monte-Carlo  method.  We  remark  that  in  dimension  d  —  3  the  speed-up  of  the  methods 
becomes  really  evident  for  large  values  of  N.  For  example,  for  N  =  64  and  M  =  4  the 
fast  method  is  more  than  14  times  faster  then  the  direct  algorithm. 


3  Asymptotic  preserving  splitting  methods 

In  this  section  we  will  consider  the  problem  of  numerically  solving  the  full  nonhomoge- 
neous  equation.  The  approximation  of  the  velocity  space  has  been  extensively  discussed 
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Figure  10:  Comparison  between  the  fast  and  classical  spectral  methods  and  the  Monte- 
Carlo  methods  for  the  temperature  components  relaxation. 

in  Section  2,  thus  here  we  will  focus  on  the  space  and  time  discretizations.  Besides  the 
solution  of  the  transport  part,  when  dealing  with  the  full  problem  the  main  difficulty  is 
the  construction  of  schemes  which  are  robust  in  the  fluid  limit,  or,  as  defined  in  Sec¬ 
tion  1.7,  which  are  asymptotic  preserving  (AP).  Of  course  AP  is  an  important  property 
strictly  related  to  the  stability  of  the  time  discretization  scheme  in  stiff  regimes.  Here  we 
will  treat  separately  methods  based  on  operator  splitting  from  other  time  discretizations 
which  avoid  operator  splitting.  First  we  start  from  operator  splitting  methods  which  have 
been  introduced  in  Section  1.6.  In  the  next  paragraphs  we  will  focus  on  the  space-time 
approximation  of  the  transport  step  and  on  the  time  discretization  of  the  collision  step. 

3.1  Transport  step 

The  solution  of  the  transport  equation  can  be  obtained  in  a  variety  of  ways  accordingly 
to  the  particular  application  considered  (see  for  example  LeVeque  (1992);  Filbet  et  al. 
(2001);  Roe  and  Sidilkover  (1992)  and  the  references  therein).  A  review  of  such  a  broad 
field  is  above  the  scope  of  this  paper.  Here  we  give  the  details  of  the  methods  developed 
originally  in  Filbet  et  al.  (2001)  which  has  several  nice  properties,  among  which  to  allow 
the  use  of  large  time  steps  even  for  large  velocities. 

3.1.1  Positive  and  Conservative  schemes 

Let  us  consider  the  transport  equation  written  as 

dtf  +  Vx(vf)  =  0,  V(t,x)eM+xlRd.  (51) 

Then,  the  solution  of  the  transport  equation  at  time  tn+l  reads 

f(tn+\x)  =  f{tn,x-vAt),  \/xeJRd. 
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Figure  11:  Time  evolution  of  the  kinetic  entropy  H  and  high-order  moments  AI4,  M.&  and 
M.%  of  f(t,v )  for  the  fast  and  classical  spectral  methods,  and  the  Monte-Carlo  methods. 
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3.1  Transport  step 


For  simplicity,  let  us  restrict  ourselves  to  a  one  dimensional  problem  and  introduce  a  finite 
set  of  mesh  points  {xi+i /2}?.e/  on  the  computational  domain.  We  will  use  the  notations 
Ax  =  xi+i/2  —  Xi_  1/2,  Ci  =  [xj_ i/2,Xj+i/2]  and  Xj  the  center  of  Q.  Assume  the  values  of 
the  distribution  function  are  known  at  time  tn  =  n  At  on  cells  Ci,  we  compute  the  new 
values  at  time  tn+1  by  integration  of  the  distribution  function  on  each  sub-interval.  Thus, 
using  the  explicit  expression  of  the  solution,  we  have 

)—v  At 


Ei  + 1/2 


f(tn+1,  x)dx  = 


Ei+ 1/2" 


Ei- 1/2 


Ei- 1/2 


—  D  At 


f(tn,x)dx , 


then,  setting 


<W(0  = 


Ei+ 1/2 


Ei+l/2 


—  V  At 


f(tn,x)dx , 


we  obtain  the  conservative  form 

fxi  +  l/2 


f{tn+1,  x)dx  = 


ci+ 1/2 


f(tn,x)dx  +  1/2 (t”)  -  *i+i/2(<”). 


(52) 


Ei- 1/2 


Ei- 1/2 


The  evaluation  of  the  average  of  the  solution  over  [xj_i/2,  Xi+1/2]  allows  to  ignore  fine 
details  of  the  exact  solution  which  may  be  costly  to  compute.  The  main  step  is  now  to 
choose  an  efficient  method  to  reconstruct  the  distribution  function  from  the  cell  average 
on  each  cell  C'j.  We  will  consider  a  reconstruction  via  primitive  function  preserving  pos¬ 
itivity  and  maximum  values  of  /  (Filbet  et  ah,  2001).  Let  F{tn,x)  be  a  primitive  of  the 
distribution  function  f(tn,x),  if  we  denote  by 

1  rxi+ 1/2 

/”=A x  I  , 

then  F(tn,x;+i/2)  -  F(tn,x^  1/2)  =  Ax  f?  and 

i 

F(tn,xi+ 1/2)  =  A  x^f, 


n  =  w?. 


k= 0 


First  we  construct  an  approximation  of  the  primitive  on  the  small  interval  [xj_i/2,  Xj+1/2] 
using  the  stencil  {xj_3/2,  Xj_i/2,  xm/2,  xi+3/2} 

Fh(r,  x)  =  <_!  +  0  -  Xi-1/2 )/,”  +  -  ^—1/2)0  -  ^+i/2)[/A-i  -  fi] 

{x  -  Xi-1/2)(x  -  Xi+1/2)(x  -  Xi+3/2)[/t"  1  -2  /f  +  /”_  J, 


6Ax2 


where  we  use  the  relation  wf  —  j  =  Ax  /”.  Thus,  by  differentiating  Fh(x),  we  ob¬ 
tain  a  third  order  accurate  approximation  of  the  distribution  function  on  the  interval 
[^1-1/2)  ^i+1/2] 

r)  fi 

A(r,x)  =  ^(r,x)  =  /r+ 

2  (x  -  Xi)(x  -  Xi_3/2)  +  (x  -  Xi_i/2)(x  -  Xi+i/2)  (/,"  1  -  fi) 

2  (x  -  Xj)  (x  -  Xj_|_3/2)  +  (x  -  Xj_i/2)  (x  -  xm/2)  (/"  -  1). 


6  Ax2  L 
1 

6  Ax2 
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Unfortunately,  this  approximation  does  not  preserve  positivity  of  the  distribution  function 
/.  Then,  in  order  to  satisfy  a  maximum  principle  and  to  avoid  spurious  oscillations  we 
introduce  slope  correctors 


fh{tn,x)=f?+ 


(53) 


+ 


£■ 


with 


6  Ax2 

0 

6  Aa;2  L 


et  = 


2{x-  Xi){x  -  Xi- 3/2)  +  (x  -  Xi-1/2) (x  -  Xi+1/2)  (/,"  1  -  /") 
2(x-  xt)(x  -  ^i+3/2)  +  (a:  -  xt-  1/2)(x  -  xi+1/2)  (/"  -  f?_  1), 


mm(l,2/,V(/&1~/") 


if  f?±1  -  f T  >  0, 


miu(i,  -2  m  -  /r)/(/;±i  -  /d)  ^  /"±i  -  /;*  <  0. 


(54) 


where  /oo  =  max{/"}  is  a  local  maximum. 

jei  J 

The  theoretical  properties  of  this  reconstruction  can  be  summarized  by  the  follow- 
ing(Filbet  et  al.,  2001) 

Proposition  2  The  approximation  of  the  distribution  function  fh(x),  defined  by  (53)- 
(54),  satisfies 

•  The  conservation  of  the  average:  for  all  i  e  /,  f* fh(x)dx  =  A x  fi. 

•  The  maximum  principle:  for  all  x  G  (xmin,xmax),  0  <  fh(x)  <  /oo- 

Moreover,  if  we  assume  that  the  Total  Variation  of  the  distribution  function  f(x)  is 
bounded,  then  we  obtain  the  global  estimate 


\fh(x)-fh(x)\dx<4TV(f)Ax, 


where  fh  denotes  the  third  order  approximation  of  f  without  slope  corrector. 

Remark  6  If  the  solution  is  smooth,  we  can  check  numerically  that  the  scheme  is  third 
order.  In  several  dimensions  we  can  perform  reconstruction  dimension  by  dimension  using 
tensor  product. 

3.2  Time  discretization  of  the  collision  step 

Since  we  aim  at  developing  operator  splitting  AP  schemes,  the  most  natural  choice  would 
be  to  use  implicit  solvers  applied  to  the  collision  step  (32).  Unfortunately  the  use  of 
fully  implicit  schemes  for  the  full  Boltzmann  collision  integral  is  impracticable  due  to 
the  prohibitive  computational  cost  required  by  the  solution  of  the  very  large  non-linear 
algebraic  system  originated  by  the  five  fold  integral  appearing  in  Q(f,f).  Exponential 
methods  represent  a  possible  way  to  overcome  these  difficulties. 
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3.2.1  Problem  reformulation 

First  we  rewrite  the  homogeneous  equation  (32)  in  the  form 

dtf  =  (55) 

where  P(f,  f )  =  Q(f ,  f)  +  gf  and  g  >  0  is  a  constant  such  that  P(f,  f )  >  0.  Typically 
g  is  an  estimate  of  the  largest  spectrum  of  the  loss  part  of  the  collision  operator. 

By  construction  we  have  the  following 

-[  P{f,f){v)<l>(v)dv=  [  f(v)<f>(v)dv,  cp(v)  =  l,v,v2.  (56) 

d  J r3  Jr3 

Thus  P(f,  /)/ yU  is  a  density  function  and  we  can  consider  the  following  decomposition 

P(f,f)/»  =  M  +  g,  (57) 

where  M  is  the  Maxwellian  with  the  same  moments  of  /. 

The  function  g  represents  the  non  equilibrium  part  of  the  distribution  function  and 
from  the  definition  above  it  follows  that  g  is  in  general  non  positive.  Moreover  since 
P(f,  f)/[i  and  M  have  the  same  moments  we  have 


g{y)4>{y)  dv  =  0,  c j){y )  =  l,v,v2 


(58) 


The  homogeneous  equation  can  be  written  in  the  form 

dtf  =  -g+  —(M  -f)  =  ^( -  m)  +^(M-  /).  (59) 

The  above  system  is  equivalent  to  the  penalization  method  for  the  collision  operator 
recently  introduced  in  Filbet  and  Jin  (2010).  Note  that  even  if  M  is  nonlinear  in  /, 
thanks  to  the  conservation  properties  (19),  it  remains  constant  during  the  relaxation 
process.  The  main  feature  of  such  formulation  is  that  on  the  right  hand  side  we  have 
a  stiff  dissipative  linear  part  g(M  —  f)/e  which  characterizes  the  asymptotic  behavior 
of  /  and  a  stiff  non  dissipative  non  linear  part  (P(f,f)/g  —  M)/e  which  describes  the 
deviations  of  P(f,  f)/g  from  M,  or  equivalently  the  deviations  of  the  Boltzmann  operator 
from  a  BGK-like  relaxation  term. 

The  formulation  above  permits  to  use  exponential  integrators  where  the  schemes  take 
advantage  of  the  exact  solution  of  the  linear  part  (Dimarco  and  Pareschi,  2010a).  Expo¬ 
nential  methods  for  kinetic  equations  were  first  proposed  in  Gabetta  et  al.  (1997). 


3.2.2  Explicit  exponential  Runge-Kutta  schemes 

In  order  to  derive  the  methods  it  is  useful  to  rewrite  (59)  as 


d (/  -  M)e^/£ 
dt 


£ 


(60) 


The  above  form  is  readily  obtained  if  one  multiplies  (59)  by  the  integrating  factor 
exp (yuf/e)  and  takes  into  account  the  fact  that  M  does  not  depend  of  time.  A  class  of 
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explicit  exponential  Runge-Kutta  schemes  is  then  obtained  by  direct  application  of  an 
explicit  Runge-Kutta  method  to  (60).  More  in  general  we  can  consider  the  family  of 
methods  characterized  by 


rc*MA  t/efn  +  pAt  AijiliAt/e)  ( _  Mn 
£  3=1  '  11 


+  (1  -  e~Ci^t/£)  Mn,  *  = 


/ 


n+l  _ 


V 


=  e 


-vAt/efn  +  W^pAt/e)  ( _  Mr 

e  tr  v  r 


+  (1  -  e~^t,£)  Mn, 


(61) 


(62) 


where  At  is  the  time  step,  fn  =  f(tn),  Mn  =  M(tn),  Ci  >  0,  and  the  coefficients  Aij  and 
the  weights  IT)  are  such  that 


^P'(O)  —  aiji  Wi(0)=Wi,  i,j  = 


with  coefficients  alJ  and  weights  Wi  given  by  a  standard  explicit  Runge-Kutta  method 
called  the  underlying  method.  Various  schemes  come  from  the  different  choices  of  the 
underlying  method.  The  most  popular  approaches  are  the  integrating  factor  (IF)  and  the 
exponential  time  differencing  (ETD)  methods  (Maset  and  Zennaro,  2009).  Since  Mn  does 
not  depend  on  time  during  the  collision  process  in  the  sequel  we  will  omit  the  index  n. 

For  the  so-called  Integrating  Factor  methods,  which  correspond  to  a  direct  application 
of  the  underlying  method  to  (60),  we  have 


Aij{ A)  =  a^e  (Ci  Cj)x,  i,  j  =  1, . .  .,i/,  i  >  j 
Wi(X)  =  Wie~(1~Ci)x,  i  =  1, is, 


(63) 


with  A  =  fiAt/e. 

The  hrst  order  IF  scheme  reads 


r+1  = 


gA  t 


=  e 


f 


uAt 
+  - — 


p(r,fri 

fi 


M  +  1 


_  fiAt 


M, 


(64) 


which  is  based  on  explicit  Euler.  For  such  methods  the  order  of  accuracy  is  the  same  as 
the  order  of  the  underlying  method. 

The  Exponential  Time  Differencing  methods  are  strictly  connected  with  the  integral 
representation  of  (60).  In  the  general  case  the  coefficients  for  ETD  methods  have  the  form 

Aij{  A)  =  [  e{1~s)CiXpij(s)  ds ,  i,j  =  i>j 

Jo 

HA:  (A)  =  [  e{1~s)xpi(s)  ds,  i  =  l,...,zy 

Jo 

where  pt  and  p^  are  suitable  polynomials. 

The  standard  hrst  order  ETD  method  based  on  explicit  Euler  in  our  case  gives 


/ 


n+l 


_  aM.  rn  pAt 
=  e  -  fn  +  —cp 
£ 


nrvn 

p 


(65) 


where  <p(z)  =  (1  —  e  z)/z. 
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3.2.3  Time  Relaxed  methods 


A  class  of  exponential  methods  for  kinetic  equations,  the  so-called  time  relaxed  (TR) 
methods,  has  been  introduced  in  Gabetta  et  al.  (1997)  as  a  combination  of  an  exponen¬ 
tial  expansion  (or  Wild  sum)  together  with  a  suitable  Maxwellian  truncation.  In  this 
paragraph  we  show  that  these  schemes  included  already  decomposition  (59)  and  can  be 
derived  directly  from  a  suitable  Taylor  expansion  of  (60). 

To  show  this,  let  us  first  introduce  the  change  of  variables 


r  =  1  —  exp  (— fit/e), 


which,  using  the  bilinearity  of  P(f,  /),  gives  the  equation 


=  (P(U)-HM)  A  ,2-  (66) 

The  application  of  an  explicit  Runge-Kutta  scheme  to  (66)  with  time  step  At  —  1  — 
exp(— /iA t/e)  leads  to  a  class  of  ETD  methods.  For  example  the  first  order  scheme  based 
on  explicit  Euler  in  the  original  variables  yields 


d 

dr 


(/  -  M) 


1  —  r 


/»+'  =  e-“  / 


fiAt 


1 1  At 

+  - — <pi 


paa  fp(r,r ) 


M  +  (1 


fxAt 

e"V)M, 


(67) 


where  <fk(z)  =  e~z(l  —  e~z)k/z,  k  =  1,  2, _ 

Note  that  such  scheme  coincides  with  the  first  order  exponential  time  relaxed  method 
derived  in  Gabetta  et  al.  (1997)  and  differs  from  the  standard  ETD  method  based  on 
explicit  Euler.  Higher  order  ETD  methods  can  be  derived  as  well  simply  applying  higher 
order  explicit  Runge-Kutta  methods  to  (66).  Although  interesting,  here  we  do  not  explore 
further  this  class  of  schemes. 

Now  let  us  consider  a  different  approach  by  taking  the  Taylor  expansion  of  ( f—M ) /(l  — 
r)  around  r  =  0  in  (66).  This  leads  to 


(/-M)/(l-r)  =  (/0  —  M)  +  r 


P(foJo) 


M 


T 

~2 


P(P(foJo)Jo)  +  P(fo,P(foJo)) 

/i2 


-2  M 


+  0(r3) 


where  we  have  used  the  bilinearity  of  the  operator  P(f,  /). 

If  we  compute  all  the  terms  in  the  expansion  and  use  recursively  the  bilinearity  of 
P(f,  f )  we  can  state  the  following 


Proposition  3  The  solution  to  problem  (59)  or  equivalently  (60)  or  (66)  can  be  repre¬ 
sented  as 

OO 

f(v,  t)  =  (1  -  T)f0(v)  +  (1  rk(fk(,v)  -  M{y))  +  tM (v),  (68) 

k=  1 

where  the  functions  fk  are  given  by  the  recurrence  formula 

k 

fk+i(v)  =  j^—^2-P(fh,fk-h)(v),  A;  =  0,1,....  (69) 

+  1  h= 0  ^ 
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By  truncating  expansion  (68)  at  the  order  m,  and  reverting  to  the  old  variables,  we 
recover  exactly  the  time  relaxed  schemes  presented  in  Gabetta  et  al.  (1997) 

m 


fn+ 1  =  e-»M/£fn  +  e~^At/e  V(1  -  e-»At/£)k(fZ  -  M)  +  (1  -  e-^/£)M, 


k= 1 


(70) 


which,  using  the  fact  that 

m 

^  t/e  Eh-  ^  [iAt/e^m+1 

k= 0 

can  be  rewritten  in  the  usual  form  emphasizing  their  convexity  properties 

m 

f  n+l  =  e-nAt/e  -  e-»A t/£)kfk  +  (1  -  e-^At/e)m+1  M.  (71) 

k= 0 

A  remarkable  feature  of  these  methods  is  that  the  functions  fk{v)  are  density  functions 
with  the  same  moments  of  the  initial  data.  Such  property,  together  with  unconditional 
nonnegativity  and  convexity  of  the  weights,  is  enough  to  guarantee  asymptotic  preser¬ 
vation  of  the  schemes  as  well  as  nonnegativity  and  entropic  stability  (see  Gabetta  et  al. 
(1997)  for  details). 


3.2.4  Main  properties 

In  this  section  we  will  report  the  main  properties  for  an  IF  exponential  scheme  in  the 
form  (61)-(62).  We  refer  to  Dimarco  and  Pareschi  (2010a)  for  further  details  an  results 
concerning  general  exponential  schemes. 

Now  let  us  denote  by  fn  and  gn  the  corresponding  solutions  obtained  with  an  explicit 
exponential  Runge-Kutta  method.  Let  us  define 


u-l 

R{ A)  =  e~A  +  ^k+1w{X)TA(X)kE{X)er  (72) 

k= 0 

where  A  =  pAt/e,  A(A)  the  v  x  v  matrix  of  elements  |Ay(A)|,  w( A)  the  v  x  1  vector  of 
elements  |W7(A)|,  e  the  v  x  1  unit  vector  and  E( A)  =  diag(e_ClA, . . .  ,e_Cl/A). 

We  can  state  (Dimarco  and  Pareschi,  2010a) 

Theorem  3  If  an  explicit  exponential  Runge-Kutta  method  in  the  form  (61)-(62)  satisfies 

lim  R( A)  =  0,  (73) 

A — ^oo 

with  R( A)  given  by  (72)  then  it  is  asymptotic  preserving. 

Note  that  for  an  IF  method  we  have 

|AAA)|  <  K-|e-(c>-c^A,  |Wj(A)|  <  \Wi\e-^\ 
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thus  we  require 

0  =  ci  <  c2  . . .  <  cv  <  1, 

in  order  for  the  above  quantities  to  be  bounded  independently  of  A. 
Moreover  for  nonnegative  coefficients  and  weights  we  get 

R{ A)  =  e~x  ( 1  +  A k+1wTAke 
V  k= 0 


(74) 


(75) 


and  condition  (73)  is  always  satisfied. 

It  can  be  proved  that  if  the  underlying  Runge-Kutta  method  is  a  //-stage  explicit 
Runge-Kutta  method  of  order  is,  with  nonnegative  coefficients  and  weights  satisfying  (74), 
then  the  scheme  is  unconditionally  stable  and  contractive.  As  pointed  out  in  Maset  and 
Zennaro  (2009)  examples  of  such  methods  are  well-known  up  to  is  —  4  and  the  classical 
RK  method  of  order  four  is  the  sole  method  with  four  stages. 

For  practical  applications  it  may  be  convenient  to  require  that  as  A  — »  oo  the  numer¬ 
ical  solution  f  n+1  and  each  level  of  the  IF  method  are  projected  towards  the  local 
Maxwellian.  It  is  straightforward  to  verify  that  this  stronger  AP  property  is  satisfied  if 
we  replace  condition  (74)  by 


0  =  Ci  <  c2  <  . . .  <  cv  <  1.  (76) 

We  conclude  this  section  with  a  results  concerning  an  important  convexity  property 
of  IF  schemes.  We  can  state  (see  Dimarco  and  Pareschi  (2010a)) 

Proposition  4  An  explicit  IF  method  is  unconditionally  positive  and  entropic  if  the  un¬ 
derlying  Runge-Kutta  method  has  nonnegative  coefficients  and  weights  satisfying  By  Tay¬ 
lor  expansion  we  obtain  conditions 

i—  1 

<  k(^v  k  =  0,1,2,...,  1  =  1,. ..,is  (77) 

3= 1 

v  X 

k  =  0,1,2,...,  (78) 


Note  that  the  above  conditions  on  the  choice  of  the  underlying  method  are  quite  restrictive 
since  we  are  not  using  the  bilinearity  of  P(f,  f)  which  would  lead  to  weaker  constraints 
on  a.ij  and  wt .  Examples  of  methods  that  satisfy  convexity  are  the  second  order  modified 
Euler  method 

jpW  —  jn 

F®  =  e~A/2fn  +  ^e~A/2  (F(F(1)^F(1))  _  M^j  +  (x  _  e-A/2)  (79) 

fU+1  =  e-A jn  +  Xe-X/2  ^P(F<-2\FW)  _  ^  +  _  e-A^  M 
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and  the  third  order  Henn  method 

jpW  —  fn, 


f 


pi2) 

_p(3) 

n+1 


=  e 


=  e 


■A/3/”  +  le“V3  All  _  A/  j  +  (1  -  e“V3)  M, 

-2A/S f„  +  ^e-A/3  ^P(f(2\f(2))  _  v/j  +  (J  _  e-2V3)  A/, 


(80) 


-  M 


e-y  +  W  (^'^(11) 

bnt  not  the  classical  fourth  order  Runge-Kutta  scheme. 

3.2.5  Implementation  aspects 

An  essential  aspect  in  the  reformulation  of  the  problem  given  by  (75)  is  the  choice  of 
the  value  of  the  constant  /i  used  in  estimating  the  spectrum  of  the  collision  operator.  Of 
course  such  constant  can  be  chosen  at  each  time  step  in  order  to  improve  our  estimate.  In 
the  sequel  we  show  different  choices  in  the  case  of  variable  hard  spheres.  We  set  C1  =  1 
for  simplicity. 

The  choice  of  an  upper  bound  for  the  loss  part  of  the  collision  term  leads  to  take 
H  —  Up  where 


Up  =  sup 

V 


\v  —  v*\ 7/(v*)  du* 


(81) 


Positivity  is  guaranteed  since  it  implies  clearly  P(f,  f)  >  0.  From  a  practical  viewpoint 
computation  of  (81)  can  be  done  at  0(N  log  N)  for  a  deterministic  method  based  on  N 
parameters  for  representing  f(v)  on  a  regular  mesh.  This  can  be  done  using  the  FFT 
algorithm  thanks  to  the  convolution  structure  of  the  loss  term  in  (81). 

However,  such  positivity  constraint  on  P(f,  f )  typically  leads  to  overestimates  of  the 
true  spectrum  of  the  collision  operator,  especially  in  Monte  Carlo  simulations.  A  better 
estimate  of  /i  would  be  to  use  the  average  collision  frequency 


\v  -  V*| 7/(n)/(n*)  ch;*  dv. 


(82) 


Finally,  as  suggested  in  Filbet  and  Jin  (2010),  //  can  be  chosen  as  an  estimate  of  the 
spectral  radius  of  the  linearized  operator  Q  around  the  Maxwellian  M.  In  fact 

Q(f,f)  »  Q(M,M)  +  VQ(M,M)(M  -  /)  =  VQ(M)(M  -  /), 

where  VQ(M,  M)  is  the  Frechet  derivative  of  Q  evaluated  at  M.  For  example  one  can 
take 

I  Q(f,f) 


\  J  ")  J  J  /  o  o  \ 

=  s“p  /Vm  •  <83> 

The  choices  /i  —  fJLa  or  /i  —  /j,s,  although  more  accurate,  pose  the  question  of  stability 
of  the  resulting  scheme  since  they  do  not  guarantee  P(f,  f)  >  0.  We  refer  to  Dimarco 
and  Pareschi  (2010a)  for  more  results  in  this  direction. 


RTO-EN-AVT-194 


8-45 


3  AP  SPLITTING  METHODS 


3.3  Numerical  tests 


Remark  7  Here  we  have  assumed  //  constant  during  the  time  stepping.  In  principle  one 
can  take  /i  —  /i(t)  and  rewrite  the  exponential  methods  for  a  time  dependent  //  from 


d(f-  M)eefo»(s)ds  1  i  ft  ,  ,  , 

ijt  £ 


(84) 


and  then  recompute  pp  at  each  time  step  or  stage  of  the  Runge-Kutta  method. 


3.3  Numerical  tests 

In  this  section  we  report  several  test  cases.  First  for  homogeneous  equations  where  we 
illustrate  the  AP  feature  of  the  exponential  schemes  and  then  for  nonhomogeneous  prob¬ 
lems  in  different  regimes. 

3.3.1  Homogeneous  problems 
A  simple  test  case 

In  the  first  test  case  we  consider  the  simplified  situation  of  the  Kac  equation  (17)  with 
initial  data  given  by 

f(v,0)=v2e~*. 

In  this  case  we  have  an  exact  solution  given  by 

/(»,  t)  =  \  0(1  -  C(t))PW)  +  (3 C(t)  -  1  )C(()3/vI  e-°«>y  (85) 

with  C(t)  =  (3  —  2e-v^t/16)-1.  Thanks  to  its  simple  structure  the  collision  integral  can 
be  evaluated  analytically  and  no  further  approximation  is  needed  when  comparing  the 
accuracy  of  different  time  discretizations  after  only  one  time  step.  More  precisely  we 
compare  the  first  and  second  order  TR  and  IF  methods  with  a  first  and  second  order 
implicit-explicit  (IMEX)  and  diagonally  implicit  Runge-Kutta  (DIRK)  methods.  The 
IMEX  methods  have  been  applied  following  Filbet  and  Jin  (2010)  (schemes  (4.3)  and 
(4.14))  which  in  a  homogeneous  setting  correspond  to  take  the  gain  part  of  the  collision 
term  explicitly  and  the  loss  part  implicitly  since  the  Maxwellian  term  cancels.  All  schemes 
are  unconditionally  stable  and  AP  except  for  the  particular  IMEX  schemes  which  do  not 
satisfy  the  AP  property.  We  refer  to  Section  4.3  for  a  discussion  of  asymptotic  preserving 
IMEX  schemes. 

The  results  are  reported  in  Figure  12.  The  AP  feature  of  the  exponential  schemes  (both 
TR  and  IF)  and  fully  implicit  solvers  permits  to  capture  the  correct  behavior.  Quite  re¬ 
markably,  in  this  test  case,  exponential  schemes  are  more  accurate  then  the  corresponding 
fully  implicit  methods  with  the  IF  methods  being  the  most  accurate. 

Accuracy  test 

Next  we  compare  the  accuracy  for  Maxwell  molecules  and  Hard  Spheres  in  3D.  As  initial 
data  we  consider  an  equilibrium  distribution  with  temperature  T  =  6,  density  g  =  1  and 
mean  velocity  u  =  —0.5.  To  this  distribution  we  add  a  bump  on  the  right  tail  along  the 
x-axis.  The  bump  is  realized  adding  a  fraction  of  particles  in  equilibrium  state  with  mass 
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Figure  12:  Kac  equation  L\  error  of  some  first  and  second  order  time  discretization 
methods 

Qb  =  0.5  q,  mean  velocity  Ub  =  4  y/T  and  temperature  Tb  =  0.5  T  to  the  initial  Maxwellian 
distribution  function.  The  velocity  error  is  neglected  by  solving  the  collision  operator  with 
Monte  Carlo  methods  with  a  very  large  number  of  particles.  The  simulations  are  run  till 
the  equilibrium  is  approximately  reached,  which  means  t  —  0.4  in  the  case  of  Hard  Spheres 
(HS)  and  t  =  0.8  for  Maxwellian  molecules  (MM).  The  reference  solution  is  computed  by 
the  same  method  with  a  very  small  time  step. 

In  Figure  13  we  show  the  L2  error  for  the  fourth  order  moment  of  the  distribution 
function  /  for  Maxwellian  molecules.  In  Figure  14  the  L2  error  for  the  fourth  order 
moment  of  /  is  reported  in  the  case  of  hard  sphere  molecules.  Observe  that,  in  the  case 
of  Maxwellian  molecules  ji  =  1  while  in  the  case  of  Hard  Sphere  particles  /x  is  a  constant 
upper  bound  for  the  collisional  cross  section  (/x  3>  1).  This  constraint  implies  that  the  L2 
norm  of  the  error  is  larger,  for  equals  choices  of  /j,At/e,  in  the  case  of  HS  respect  to  the 
case  of  Maxwellian  molecules. 
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Figure  13:  L2  Error  for  the  Fourth  Moment  Relaxation  for  the  homogeneous  relaxation 
problem  with  Maxwellian  particles. 
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Figure  14:  L2  Error  for  the  Fourth  Moment  Relaxation  for  the  homogeneous  relaxation 
problem  with  hard  sphere  particles. 
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Adaptivity 

In  this  test  we  explore  the  possibility  to  use  an  adaptive  value  for  /j  in  time.  We  use  a 
first  order  TR  scheme  with  adaptive  time  stepping  and  /i  given  by 

f*(t)  =  3/2  max  ^L(v)  -  , 

where  fl  =  [— T/2,T/2]3.  We  use  Maxwell  molecules  with  the  spectral  method  and  642 
modes  in  2D  so  that  the  velocity  error  can  be  neglected.  We  compare  the  results  with 
a  reference  solution  obtained  with  a  4th  order  Runge-Kutta  with  a  fixed  time  step.  The 
results  are  reported  in  Table  2. 

The  time  evolution  of  //  is  given  in  Figure  15.  It  shows  that  the  method  becomes  more 
and  more  accurate  as  the  solution  approaches  the  Maxwellian  state. 

3.3.2  Non-homogeneous  problems 
Accuracy  of  the  method 

Here  we  consider  the  full  non-homogeneous  case  where  the  collision  operator  is  solved  by 
spectral  methods,  the  transport  by  the  flux-positive  schemes  and  the  operator  splitting  is 
second  order  Strang  splitting  combined  with  the  adaptive  second  order  TR  method.  We 
test  the  overall  accuracy  of  the  method  using  as  initial  condition 

fo(x,v)  =  (1  +  {3  cos(kox))  exp(—v2/2),  (x,v)  G  [0 ,L\  x  JR2, 

with  periodic  boundary  conditions.  The  error  is  computed  as 

e2h  =  max  (||  fh{t)  -  /2h(t)||i)/||/o||i, 
te(o,T) 

and  the  results  are  given  in  Table  3.  As  expected  the  second  order  accuracy  of  the  scheme 
is  observed. 


r  =  n(t)  At 

nrot 

Numerical  error  A1) 
with  a  fixed  At  =  T /riTot 

Numerical  error  e ^ 
with  At  =  r/n(t) 

0.010 

50 

0.0036 

0.002 

0.025 

20 

0.0080 

0.007 

0.050 

11 

0.0160 

0.014 

0.100 

08 

0.0300 

0.025 

0.200 

07 

0.0520 

0.045 

0.500 

05 

0.2000 

0.090 

1.000 

03 

XXXX 

0.150 

3.000 

02 

XXXX 

0.040 

5.000 

01 

XXXX 

0.006 

Table  2:  Comparison  between  fixed  A t  and  At  =  r/n(t) 
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Figure  15:  Time  evolution  of  the  value  p,(t). 


Numerical  parameters 

relative  l 1  error  norm 

£2 h/th 

nx  =  032,  nVx  =  nVy  =  08,  At  =  0.100 

£4  h  =  0.3835 

5.20 

nx  =  064,  nVx  =  nVy  =  16,  At  =  0.050 

e2h  =  0.0738 

4.35 

nx  =  128,  nVx  =  nVy  =  32,  At  =  0.025 

£h  =  0.0169 

X 

Table  3:  Convergence  results 


Stationary  shock  profile 

We  consider  a  stationary  shock  wave  problem  for  the  Boltzmann  equation  solved  on  a 
finite  domain  —  L  <  x  <  L  with  boundary  conditions  that  the  incoming  flux  of  particles 
at  x  —  ±L  is  distributed  according  to  the  Maxwellian  flux  vM±{y).  As  initial  data,  we 
take  f(x,v,  0)  =  M(p,u,T),  with 

p  =  1.0,  T  =  1.0,  M.  =  2.0,  L  >  x  >  0, 

where  At  is  the  Mach  number.  The  mean  velocity  is 

ux  =  —My/jT,  uy  =  0, 


with  7  =  5/3. 

The  values  for  p,  u  and  T  for  x  <  0  are  given  by  the  Rankine-Hugoniot  conditions 
(Whitham,  1974). 

The  profiles  are  shown  in  Figure  16  for  different  Knndsen  numbers.  As  a  reference 
solution  we  report  also  the  solution  obtained  by  Monte  Carlo  methods. 

Riemann  problem 

This  test  deals  with  the  numerical  solution  of  the  non  homogeneous  IP  x  2 D  Boltzmann 
equation  for  hard  sphere  molecules  (a  =  1).  We  have  computed  an  approximation  for 
different  Knudsen  numbers,  from  rarefied  regime  up  to  the  fluid  limit.  The  solution  in  the 
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Figure  16:  Shock  profiles  at  Mach  2:  £  =  10  1  (left)  and  £  =  0.05  (right).  From  top  to 
bottom:  density  p,  mean  velocity  u  and  temperature  T. 
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hydrodynamic  limit  is  also  compared  with  the  numerical  solution  of  Euler  system.  The 
initial  data  is  given  by 

(  (phuhTl)  =  {  1,0,1)  if  0  <  x  <  0.5, 

(  (pr,  ur,  Tr)  =  (0.125,  0,  0.25)  if  0.5  <  x  <  1, 

In  Fig.  17  we  plot  the  results  obtained  in  the  rarefied  regime  (e  =  10”2)  using  the 
Spectral-PFC  scheme  and  a  Monte  Carlo  method  (TRMC)  as  a  comparison.  The  TRMC 
method  is  used  with  100  cells  in  x  containing  100  particles  whereas  the  Spectral-PFC 
scheme  is  used  with  64  points  in  x  and  the  size  of  the  velocity  grid  is  64  x  64  points  for 
the  transport  and  the  total  number  of  modes  32  x  32.  We  observe  that  the  two  solutions 
are  in  this  case  very  comparable  even  if  small  oscillations,  due  to  the  statistical  noise, 
persist.  We  also  give  the  result  of  the  computations  close  to  the  Euler  limit  (e  =  10~4) 
using  128  space  cells  for  the  Spectral-PFC  method. 

Finally,  the  profiles  obtained  with  TRMC  and  Spectral-PFC  methods  are  reported  in 
Fig.  18.  On  the  opposite,  using  a  small  time  step  (A t  =  0.001),  an  accurate  solution  is 
obtained  by  the  Spectral-PFC  method,  which  is  much  less  diffusive  then  the  Monte  Carlo 
methods. 


S-PFC 

•""-v  TRMC 


0.2  0.4  0.6  0.8 


/  \ 

/  \ 

0.2  0.4  0.6  0.8  1 


0.2  0.4  0.6  0.8  1 


(1)  (2)  (3) 

Figure  17:  Riemann  problem  (kn  =  10~2):  evolution  of  (1)  the  density  p,  (2)  mean  velocity 
u  and  (3)  temperature  T  at  t  =  0.05,  0.15,  0.20. 
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(1)  (2)  (3) 

Figure  18:  Riemann  problem  (kn  =  10-4):  (1)  the  density  p,  (2)  mean  velocity  u  and  (3) 
temperature  T  at  t  =  0.20  obtained  by  the  central  scheme  for  Euler  equations  (up)  and  by 
Spectral-PFC  and  TRMC  methods  for  Boltzmann  equations. 


The  ghost  effect 

Consider  a  gas  between  two  plates  at  rest  in  a  finite  domain.  In  this  situation,  the 
stationary  state  at  a  uniform  pressure  (the  velocity  is  equal  to  zero  and  the  pressure  is 
constant)  is  an  obvious  solution  of  the  Navier-Stokes  equations;  the  temperature  held  is 
determined  by  the  heat  conduction  equation 

u  =  0,  T  =  C  -  Vx(T1/2V.tT)  =  0. 

On  the  other  hand,  if  we  move  the  plate  by  a  velocity  proportional  to  the  Knudsen  number, 
then  the  macroscopic  fields  (density  and  temperature  profiles)  will  be  affected  by  the  how, 
even  for  vanishing  Knudsen  number.  This  effect,  called  ’’ghost  effect”,  is  predicted  by  the 
Hilbert  expansion  of  the  Boltzmann  equation  in  terms  of  the  Knudsen  number,  and  it  is 
rather  difficult  to  capture  numerically,  since  the  flow  velocity  is  very  small  (Sone  et  ah, 
1996).  The  results  show  that  the  numerical  solution  agrees  with  one  obtained  by  the 
asymptotic  theory  and  not  with  the  one  obtained  from  the  heat  conduction  equation;  this 
result  is  a  confirmation  of  the  validity  of  the  asymptotic  theory. 

Thus  consider  two  parallel  plates,  both  with  temperature  distribution 

Tw(x)  =  1  —  0.5  cos(2  7ra;);  Wx  G  (0, 1), 
in  slow  motion  with  velocity 

uw(x)  =  (e,  0). 

We  use  the  hard  spheres  model  with  diffusive  b.c.  on  the  walls  and  periodic  in  x.  The 
cross  section  of  temperature  and  velocity  profile  are  shown  in  Figure  19  for  various  values 
of  the  Knudsen  number,  while  velocity  field  and  isothermal  lines  are  reported  in  Figure 
20,  for  Knudsen  number  e  =  0.02. 
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0  0.2  0.4  0.6  0.8  1  0  0.2  0.4  0.6  0.8  1 


Temperature  at  y  =  0.5  Mean  velocity  u  at  y  —  0.5 

Figure  19:  Ghost  effect:  temperature  and  mean  velocity  along  y  =  const  for  various 
Knudsen  numbers  e  =  0.05,0.02,0.01. 


velocity  field  u  isothermal  lines 


Figure  20:  Ghost  effect:  velocity  held,  and  isothermal  lines;  Knudsen  number  £  =  0.02. 
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4  Other  asymptotic  preserving  methods 

In  this  section  we  shall  consider  alternative  approaches  to  splitting  methods  when  dealing 
with  the  time  integration  of  kinetic  equations  with  stiff  collision  operators.  This  could 
be  an  advantage  for  the  construction  of  high-order  or  well-balanced  schemes,  i.e.  schemes 
that  preserve  the  stationary  solutions.  First  we  discuss  semilagrangian  schemes  and  then 
IMEX  schemes  for  BGK  like  models.  As  we  shall  see,  the  BGK  model  allows  a  simple 
implicit  treatment.  This  is  useful  per  se,  and  it  provides  a  building  block  for  effective 
treatment  of  the  Boltzmann  equation  near  the  fluid  dynamic  limit. 

4.1  Semilagrangian  methods  for  the  BGK  model 

Here  we  will  focus  on  implicit  semilagrangian  scheme  for  the  numerical  solution  of  the 
BGK  model  of  the  Boltzmann  equation.  We  shall  restrict  to  the  BGK  equation  in  one 
space  dimension.  More  details  on  the  method  can  be  found  in  Santagati  (2007). 

4.1.1  A  basic  first  order  scheme 

Let  us  rewrite  the  BGK  model  in  one  dimension: 

|+„|M(M[/]-/).  (86) 

The  numerical  scheme  for  the  solution  of  Eq.  (86)  is  based  on  the  characteristic  for¬ 
mulation  of  the  problem  (86), 

\(M  [/]-/), 

(87) 

v,  v  ' 

x,  /(0,  x,  v)  =  /o(5,  v)  t>  0, 

For  simplicity  we  assume  constant  time  step  At  and  uniform  grid  in  physical  and  velocity 
space,  with  mesh  spacing  Ax  and  Av  respectively,  and  denote  the  grid  points  by  tn  =  nAt, 
Xi  =  iAx,  i  =  1,  •  •  • ,  Nx,  Vj  =  j Av,  j  =  —Nv, . . . ,  Nv,  where  Nx  and  2 Nv  +  1  are  the 
numbers  of  grid  nodes  in  space  and  velocity,  respectively.  Let  /”  denote  the  approximate 
solution  of  the  problem  (87)  at  time  tn  in  each  spatial  and  velocity  node. 

We  start  by  considering  first  order  accurate  schemes. 

An  explicit  first  order  semilagrangian  scheme  could  be  constructed  by  computing  an 
approximation  /  of  f(tn+1,Xi  +  VjAt,Vj)  as 

+  v,A t,  v,)  =  /■>  +  -  /£)  (88) 

The  function  /  computed  in  this  way  at  the  new  time  step  does  not  lie  on  a  grid.  The  values 
of  /”+ 1  could  be  reconstructed  from  the  computed  values  /  by  a  suitable  interpolation 
back  on  the  grid  points  (see  Figure  21).  Let  us  denote  by  xt  —  ay  +  VjA t,  and  let  us 
assume  that  this  point  is  between  point  Xk  and  Xk+i •  Then  the  function  /^+1  can  be 


(l  = 

dt 

dx 

dt 

x(0)  = 
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t 


rn+l 


Figure  21:  Propagation  of  the  information  along  characteristics  in  the  explicit  scheme. 
The  computed  point  does  not  lie  on  the  grid,  and  some  interpolation  is  needed  in  order 
to  compute  /)”+ 1 .  In  this  case  k  =  i  +  2. 

reconstructed  by  simple  linear  interpolation  using  the  computed  value  at  points  Xj  and 

Xi- 1- 

The  Maxwellian  =  M(vj,  {p",w",T”})  is  computed  as  follow 


(89) 


This  formula  requires  the  computation  of  the  discrete  moments  of  {/,”}.  This  can  be 
done  by  using  a  numerical  approximation  of  the  integrals  computed  in  (21).  Following 
the  notation  in  Mieussens  (2000),  the  discrete  velocity  grid  may  be  denoted  by  V,  which 
is  composed  of  2 Nv  +  1  nodes,  and  the  moments  of  any  quantity  g,  <  g  >,  can  be 
approximated  by  a  quadrature  rule  on  V.  Let  <  g  >v  denote  the  approximation  of 
<  g  >,  where  V  is  the  set  of  2 Nv  +  1  indices  matching  the  velocity  grid  nodes.  By  this 
way  we  compute  the  moments  of  the  Maxwellian  at  each  grid  nodes  {x,}, 

(pi,PiUi,Ei)  =<  fi<j>{v)  >v 

As  quadrature  rule  we  use  summation  over  V  times  Av,  providing  spectral  accuracy  for 
smooth  functions  on  compact  support.  The  grid  V  is  chosen  to  include  most  of  the  mass. 
For  a  given  number  of  nodes  Nv,  an  optimal  choice  of  the  grid  is  obtained  as  a  compromise 
between  the  extension  of  the  velocity  domain  and  the  resolution  of  the  grid. 

Once  the  moments  are  computed  on  the  grid,  they  can  be  in  turn  computed  in  Xj, 
by  a  suitable  interpolation  formula,  so  that  the  Maxwellian  gets  easily  evaluated.  Details 
about  the  WENO  interpolation  can  be  found  in  Cockburn  et  al.  (1998),  or  in  Russo  and 
Santagati  (2011). 

The  scheme  (88)  can  be  used  to  perform  the  time  step.  The  scheme  is  first  order 
accurate.  A  more  sophisticated  time  integrator,  coupled  with  a  more  accurate  interpola¬ 
tion,  would  provide  greater  accuracy  in  time.  Notice  that,  because  of  the  semilagrangian 
nature  of  the  method,  there  is  no  CFL-type  stability  restriction  on  the  time  step  due  to 
convection.  However,  such  scheme  would  suffer  from  stability  restriction  on  the  time  step 
due  to  the  collision  term  when  the  collision  time  e  is  small. 

In  order  to  circumvent  the  stiffness  arising  from  when  e  is  small,  it  is  possible  to  resort 
to  an  implicit  formulation. 
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By  applying  simple  implicit  Euler  on  the  characteristic  equation  in  order  to  compute 
/T+1  one  obtains: 

/”+1  =  /« + f(M5+1  “  -C1)  (9°) 

The  quantity  /”■  «  f(tn,Xi  —  VjA t,Vj)  can  be  computed  by  suitable  reconstruction 
from  {/"};  linear  reconstruction  will  be  sufficient  for  first  order  scheme,  while  higher 
order  reconstruction,  such  as  ENO  or  WENO,  could  be  used  to  provide  higher  order  non 
oscillatory  reconstruction. 

The  equation  cannot  be  immediately  solved  for  f£+1,  because  the  Maxwellian  depends 
from  fn+1  itself.  However,  one  can  act  as  follows:  let  us  take  the  moments  of  both  sides 
of  Eq.  (90).  This  is  obtained  at  a  discrete  level  multiplying  both  sides  by  (j)j An,  where 
<f>j  =  {l,Vj,\vj\2},  and  summing  over  j.  We  denote  this  procedure  by  <  •  >y,  i.e.  for  any 
quantity  hj ,  we  define 

<  h  >v=  hjAv 

VjZV 

Then  we  have 

<  (ft'  -  f~t<h  >v=  Y  <  (M«+1  -  ft't  >v  (91) 

The  right  hand  side  is  zero,  because,  by  definition,  the  Maxwellian  at  time  tn+1  has  the 
same  moments  as  fn+1.  As  a  consequence,  the  moments  of  fn+1  can  be  computed  as 
moments  of  ftJ .  More  explicitly,  one  has 


p't 

= 

vjeV 

<+1 

= 

1  *  Vjev 

(92) 

1  i 

=  “•)2A‘- 
1  1  Vjev 

where  d  denotes  the  dimension  of  velocity  space  (in  our  case  d  —  1). 

Once  the  density,  mean  velocity  and  temperature  are  computed,  then  the  Maxwellian 
at  the  new  time  step  can  be  explicitly  computed: 

Once  the  Maxwellian  is  known,  the  distribution  function  can  be  explicitly  computed: 


f  f n.  +  AtMn+1 

rn+ 1  _  ^ 

Jij  e  +  At 


(93) 


Notice  that  since  we  use  a  semilagrangian  scheme,  the  geometrical  CFL  condition  is 
always  satisfied,  since  the  value  /”■  is  computed  by  interpolation  from  points  that  surround 
points  Xij  =  Xi  —  VjAt.  For  example,  in  the  case  illustrated  in  Figure  22,  xtj  e  [xk,Xk+ 1], 
with  k  =  i  —  2,  and  for  the  first  order  scheme  the  value  of  /()  can  be  obtained  by  linear 
interpolation  from  the  values  and  fg+1  •. 
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Figure  22:  Propagation  of  the  information  along  characteristics  in  the  implicit  scheme. 
The  foot  of  the  characteristics  does  not  lie  on  the  grid,  and  some  interpolation  is  needed 
in  order  to  compute 


Remark  8  Since  the  moments  are  computed  by  a  quadrature  formula,  it  is  not  properly 
true  that,  in  the  discrete  formulation,  M[f ]  and  f  have  the  same  moments.  To  get  an 
insight  on  this  aspect  see  Mieussens  (2000).  In  that  paper  the  author  introduces  the  notion 
of  a  discrete  Maxwellian,  which  is  more  consistent  with  the  discrete  formulation  of  the 
problem.  The  discrete  BGK  model  obtained  using  such  Maxwellian  is  conservative  and 
entropic.  By  enough  large  number  of  grid  points  in  velocity,  the  continuous  and  discrete 
Maxwellians  give  comparable  results.  However,  for  coarse  discretization  in  velocity,  the 
discrete  Maxwellian  introduced  in  Mieussens  (2000)  produces  better  results. 


In  next  section  we  shall  show  how  to  generalize  the  procedure  and  construct  higher 
order  schemes. 

4.1.2  High  order  time  discretization 

System  (87)  is  a  typical  ordinary  differential  equation  with  relaxation,  to  be  solved  in  the 
characteristics  framework.  Relaxation  time  lies  in  a  very  wide  range.  It  typically  extends 
from  order  one  to  very  small  values  compared  to  the  time  scale  of  the  problem.  For 
this  reason,  we  treat  the  relaxation  operator  by  L-stable  diagonally  implict  Runge  Kutta 
(DIRK)  schemes  (Hairer  and  Wanner,  1996;  Pareschi  and  Russo,  2000a,  2005),  which 
provides  enough  stability  to  ensure  the  AP  property  (see  Sec.  1.7).  DIRK  schemes  are 
defined  by  the  triangular  zzx  v  matrix,  A  =  (a^),  and  the  coefficient  vectors,  c  =  (1, ...,  cv)T 
and  b  =  (1, ...,  v)T ,  which  are  derived  by  imposing  accuracy  and  stability  constraints.  They 
characterize  completely  a  DIRK  scheme,  which  can  be  rappresented  by  the  Butcher’s 
tableaux 


c 

A 

wT 

The  internal  stages  are  practically  evaluated  by  a  sequence  of  elementary  implicit  Euler 
steps,  because,  at  each  stage,  only  the  last  stage  value  is  unknown,  due  to  the  triangular 
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4.1  Semilagrangian  methods 


structure  of  matrix  A.  Scheme  (90)  corresponds  to  implicit  Euler  scheme.  It  will  be 
denoted  Ml.  The  DIRK  methods  considered  in  this  work  are 


a 

a 

M2  =  1 

1  —  a 

a 

1  —  a 

a 

M3 


1/2 

7 

(l  +  7)/2 

(l-7)/2  7 

1 

1  —  5  —  7  <5  7 

1  —  5  —  7  <5  7 

which  are  a  second  and  a  third  order  L-implicit  schemes  Pareschi  and  Russo  (2000a).  The 
coefficients  are 


a  = 


7  =  0.4358665215,  5 


-0.644373171. 


Some  DIRK  schemes  have  the  property  that  the  last  row  of  matrix  A  equals  the  vector 
of  the  weights  b.  Such  schemes  are  called  stiffly  accurate.  This  property  is  related  to  the 
L-stability  of  the  scheme.  For  more  details  see,  for  example,  the  classical  book  by  Hairer 
and  Wanner  (1996). 

Here  we  apply  the  Runge-Kutta  scheme  along  the  characteristics.  The  numerical 
solution  is  obtained  as 

V 

/5+1  =  ;iJ  +  A(5>4')  (94) 

1=1 

The  quantity 

7?  =  )(v[if  >]  -  if) 

are  the  so  called  Runge-Kutta  fluxes,  and  have  to  be  evaluated  along  the  characteristics, 
and  depend  on  the  so  called  stage  values  F-? .  In  a  standard  DIRK  methods,  the  Ath 

stage  value,  say  F?\  is  evaluated  by  solving  an  implicit  equation  involving  only  F-? ,  since 
the  previous  stage  values  have  already  been  computed,  due  to  the  triangular  structure  of 
the  matrix  A.  In  our  case,  however,  this  is  not  so  easy,  because  if  the  point  corresponding 
to  stage  l  along  the  characteristics  is  not  a  grid  point,  it  is  not  possible  to  compute  the 
moments  of  the  Maxwellian  at  that  point  in  space-time.  For  this  reason,  one  has  to  resort 
to  a  triangular  scheme,  which  is  illustrated,  for  a  two  stage  scheme,  with  the  help  of  Figure 
23. 

Two  kinds  of  stage  values  will  be  needed:  the  stage  values  along  the  characteristics, 
denoted  by  F?\  which  are  needed  for  the  computation  of  the  numerical  solution,  and  the 

stage  values  on  the  grid,  denoted  by  Ff-  ,  which  can  be  computed  implicitly. 

The  two-stage  stiffly-accurate  scheme  M2  works  as  follows. 

First  we  compute  P[?  by 


eft?  +  A  tM, 


(i) 

ij 


e  +  At 


(95) 


Here  the  Maxwellian  M-?  =  M[P[?]  can  be  computed  by  computing  the  moments  of 
ft?  =  ffixi  —  avjAt)  by  suitable  space  reconstruction  at  time  tn. 

Once  the  implicit  step  is  solved,  the  value  of  the  Runge-Kutta  flux  k[?  can  be  easily 
computed  at  the  grid  points  ( Xi,tn  +  a  At).  Then  the  fluxes  are  computed  by  high  or¬ 
der  interpolation  on  the  intermediate  nodes  along  the  characteristics  (marked  by  a  blue 
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Figure  23:  Construction  of  two-stage  stiffly-accurate  DIRK.  The  first  stage  is  computed 
on  the  grid  points  {tn  +  aAt,xf)  (empty  red  circle)  by  an  implicit  step.  The  RK  fluxes  are 
computed  on  the  same  points,  and  then  interpolated  to  points  (tn  +  aAt,  Xj  —  (1  —  a)vj At), 
denoted  by  a  small  blue  square.  Finally,  the  second  stage  value  along  the  characteristics 
is  computed  on  the  grid.  Such  stage  value  corresponds  to  the  numerical  solution  implicit 
scheme 


square).  Once  such  values,  K-j\  are  computed,  then  the  value  of  the  numerical  solution 
can  be  computed  by 

=  %  +  M,)  (*) 

Notice  that  in  this  fl+1  =  F^\  because  the  scheme  is  stiffly  accurate,  i.e.  <221  =  &i  and 
0-22  = 

The  generalization  of  the  procedure  to  high  order  schemes  is  straightforward,  and  can 
be  found,  for  example  in  Russo  and  Santagati  (2011). 


Remark  9  In  practice  the  Runge  Kutta  fluxes  can  be  computed  from  the  internal  stages. 
For  example 


At 

£ 


Hence  the  scheme  can  be  used  in  the 
amplitude. 


p(l)  _  £71 

rij  Jij 
an 

limit  £  — >  0,  with  no  constraint  on  the  time  step 


-  FF1) 
11  ) 


4.2  Numerical  tests 

These  tests  are  aimed  to  verify  the  accuracy  (test  1)  and  the  shock  capturing  properties 
(test  2)  of  the  schemes. 

4.2.1  Regular  velocity  perturbation 

This  test  has  been  proposed  in  Pieraccini  and  Puppo  (2007).  The  solution  is  smooth,  and 
the  accuracy  can  be  tested.  Initial  velocity  profile  is  given  by 

u0(x)  =  —  (exp  (— (ax  —  l)2)  —  2 exp  (—  (ax  +  3)2))  ,  x  G  [—1, 1] 
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where  o  is  a  positive  constant  parameter.  Initial  density  and  temperature  profiles  are 
uniform,  with  constant  value,  p  —  1  and  T  =  1  respectively.  The  initial  condition  for 
the  distribution  function  is  the  Maxwellian,  computed  by  given  macroscopic  fields.  The 
boundary  conditions  are  imposed  assuming  that  beyond  the  computational  domain  the 
distribution  function  is  a  Maxwellian  with  prescribed  moments.  Two  regimes  (rarefied  and 
fluid)  have  been  investigated,  corresponding  to  different  Knundsen  numbers,  £  =  1CT2  and 
£  =  10-6.  The  final  time,  for  both  cases  was  0.04,  large  enough  to  reach  thermodynamic 
equilibrium.  Accuracy  and  conservation  tests  have  been  performed  at  the  final  time. 
The  errors  has  been  computed  using  a  reference  solution,  defined  on  a  hirer  grid,  with 
Nx  =  1280  and  Nv  =  20.  The  test  case  has  been  performed  using  Nv  =  20  (as  for  the 
reference  solution), for  each  spatial  grid  nodes  number,  uniformly  spaced  in  [-10,10]. 

The  relative  errors  and  order  of  accuracy  are  shown  in  Tables  4-7,  for  the  schemes  M2 
and  M3.  Notice  that  only  a  moderate  improvements  is  obtained  using  M3  in  place  of 
M2.  The  reason  is  the  following.  The  space  reconstruction  in  both  schemes  is  WENO 
2-3,  which  guarantees  third  order  accuracy  for  smooth  functions.  Since  space  errors  are 
dominant  with  respect  to  time  errors,  then  only  a  small  improvement  is  gained  by  a  more 
accurate  time  integrator.  This  also  explains  why  for  small  number  of  grid  points  the  order 
or  accuracy  appears  between  two  and  three. 

Also,  notice  that  for  smaller  values  of  £  a  strong  degradation  of  the  accuracy  is  ob¬ 
served.  In  fact,  in  the  fluid  dynamic  regime,  only  Erst  order  accuracy  in  time  is  guaranteed 
by  this  approach. 

Conservation  errors  have  been  investigated.  Despite  the  schemes  are  not  strictly  con¬ 
servative,  conservation  properties  are  maintained  with  good  accuracy,  even  for  a  moderate 
number  of  grid  points. 

4.2.2  Riemann  problem 

This  test  allows  us  to  evaluate  the  capability  of  our  class  of  schemes  in  capturing  shocks, 
contact  discontinuities  and  the  density  profile  in  a  rarefaction.  The  macroscopic  fields  are 
initially  assigned  in  the  domain,  satisfying  the  Rankine-Hugoniot  shock  jump  conditions. 
In  particular  we  are  interested  in  the  behavior  in  the  fluid  regime.  Here  we  illustrate  the 
results  in  the  moments,  i.e.  density,  velocity  and  temperature  profiles,  for  £  =  ICE2  and 
£  =  10~6,  respectively.  As  in  test  1,  the  boundary  conditions  are  imposed  by  Maxwellians 
computed  by  prescribed  macroscopic  moments. 

For  this  test  two  values  £  are  employed  ,  £  =  ICE2  and  £  =  ICE6.  Nv  =  60  nodes 
are  used  in  the  range  [-10,10]  of  the  discrete  velocity  domain,  as  in  Pieraccini  and  Puppo 
(2007). 

As  it  appears  from  the  results,  the  scheme  is  able  to  capture  the  fluid  dynamic  limit  for 
very  small  values  of  the  relaxation  time,  where  the  evolution  of  the  moments  is  governed 
by  the  underlying  Euler  equations. 

Remark  10  The  scheme  presented  here  is  not  in  conservation  form.  It  is  possible  to 
construct  a  conservative  version  of  the  scheme,  by  using  the  non  conservative  scheme  as 
a  predictor,  and  a  postprocessing  that  acts  as  a  conservative  corrector  to  be  used  at  each 
time  step.  The  details  of  the  scheme  will  be  fowid  in  Russo  and  Santagati  (2011). 
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L2-Relative  errors 

L 

,2-Orders 

Nx 

Density 

Velocity 

Temperature 

Dens. 

Vel. 

Temp 

20 

2.54838e-03 

2.35049e-03 

4.63423e-03 

- 

- 

- 

40 

5.57339e-04 

3.64146e-04 

9.17053e-04 

2.193 

2.690 

2.337 

80 

8.41532e-05 

5.21314e-05 

1.28062e-04 

2.727 

2.804 

2.840 

160 

1.17817e-05 

8.94658e-06 

2.53109e-05 

2.836 

2.543 

2.339 

320 

1.69746e-06 

1.95126e-06 

6.47027e-06 

2.795 

2.197 

1.968 

Table  4:  Scheme  M2,  e  —  10  2,  CFL  =  4.5. 


L2-Relative  errors 

L2-Orders 

Nx 

Density 

Velocity 

Temperature 

Dens. 

Vel. 

Temp. 

20 

2.96809e-03 

3.15227e-03 

6.37101e-03 

- 

- 

- 

40 

6.57722e-04 

6.05216e-04 

1.94043e-03 

2.174 

2.381 

1.715 

80 

1.11120e-04 

1.36059e-04 

5.26168e-04 

2.565 

2.153 

1.883 

160 

2.25137e-05 

4.61239e-05 

1.59816e-04 

2.303 

1.561 

1.719 

320 

6.10643e-06 

1.63245e-05 

5.05549e-05 

1.882 

1.498 

1.660 

Table  5:  Scheme  M2,  e  =  10  6,  CFL  =  4.5. 


L2-Relative  errors 

L2-Orders 

Nx 

Density 

Velocity 

Temperature 

Dens. 

Vel. 

Temp. 

20 

2.41539e-03 

1.97185e-03 

4.36445e-03 

- 

- 

- 

40 

4.93444e-04 

2.90747e-04 

8.14397e-04 

2.291 

2.762 

2.422 

80 

7.36995e-05 

4.27296e-05 

1.14397e-04 

2.743 

2.766 

2.832 

160 

1.06248e-05 

5.91413e-06 

1.54660e-05 

2.794 

2.853 

2.887 

320 

1.55051e-06 

9.55269e-07 

2.89414e-06 

2.777 

2.630 

2.418 

Table  6:  Scheme  M3,  £  =  10  2,  CFL  =  4.5. 


L2-Relative  errors 

L2-Orders 

Nx 

Density 

Velocity 

Temperature 

Dens. 

Vel. 

Temp. 

20 

2.02024e-03 

2.72023e-03 

4.79517e-03 

- 

- 

- 

40 

5.04007e-04 

7.57946e-04 

2.06197e-03 

2.003 

1.844 

1.218 

80 

9.91878e-05 

2.63921e-04 

9.43964e-04 

2.345 

1.522 

1.127 

160 

4.10972e-05 

1.48278e-04 

5.14779e-04 

1.271 

0.932 

0.975 

320 

2.17256e-05 

7.52320e-05 

2.39514e-04 

1.120 

1.079 

1.104 

Table  7:  Scheme  M3,  e  =  10  6,  CFL  =  4.5. 


4.3  IMEX  Runge-Kutta  schemes 

The  possibility  of  implementing  a  very  efficient  implicit  solver  for  the  BGK  model,  coupled 
with  the  fact  that  the  BGK  operator  itself  is  a  crude  approximation  of  the  Boltzmann 


8-64 


RTO-EN-AVT-194 


4.3  IMEX  Runge-Kutta  schemes 


4  OTHER  AP  METHODS 


e  =le-2 

Nx 

Density 

Momentum 

Energy 

20 

3.19103e-04 

6.45517e-04 

6.44489e-04 

40 

8.68364e-05 

2.43064e-05 

1.53479e-04 

80 

1.86619e-05 

6.93736e-06 

3.38482e-05 

160 

2.21369e-06 

6.03659e-07 

3.82991e-06 

320 

2.47089e-07 

6.48153e-08 

4.23732e-07 

640 

2.75503e-08 

5.92434e-09 

4.57648e-08 

1280 

3.21756e-09 

6.82768e-10 

5.29436e-09 

£  =le-6 

Nx 

Density 

Momentum 

Energy 

20 

3.86571e-04 

8.49545e-04 

8.59503e-04 

40 

1.25170e-04 

4.22174e-05 

2.34721e-04 

80 

3.17682e-05 

1.53056e-05 

6.50680e-05 

160 

4.65888e-06 

1.86177e-06 

9.55101e-06 

320 

5.94587e-07 

2.25448e-07 

1.20946e-06 

640 

7.47884e-08 

2.69514e-08 

1.52337e-07 

1280 

9.24536e-09 

3.27055e-09 

1.87483e-08 

Table  8:  Errors  in  conservation.  Scheme  M2. 

collision  operator,  can  be  used  as  a  tool  to  construct  very  effective  solvers  for  the  full 
Boltzmann  equation.  In  this  section  we  will  consider  this  possibility  in  the  context  of 
Implicit-Explicit  Runge-Kutta  schemes. 

4.3.1  General  formulation  of  IMEX  schemes 

To  this  goal  we  introduce  the  general  formulation  of  the  IMEX  schemes  (Pareschi  and 
Russo,  2005;  Pieraccini  and  Puppo,  2007)  for  the  BGK  model.  A  general  IMEX  schemes 
applied  to  a  kinetic  equation  of  the  form 

a,f  +  vVJ=-£(M[f]~f)  (97) 

reads 

i—  1  v  . 

=  fn  -  A  •  VXF®  +  A  ai3-{M[F^}  -  F&) 

3= 1  3= 1  ^ 

V  V  1 

fn+ 1  =  fn  -  A t^UiV  ■  VXF^  +  At^^-(M[F(i)]  -  Fw). 

*=1  j=l  5 

The  matrices  A  =  (atJ ) ,  at]  =  0  for  j  >  i  and  A  =  (atj)  are  v  x  v  matrices  such  that  the 
resulting  scheme  is  explicit  in  Vxf,  and  implicit  in  M[f]  —  f.  In  general,  an  IMEX  Runge- 
Kutta  scheme,  is  characterized  by  the  above  defined  two  matrices  and  the  coefficient 
vectors  w  =  (wi, ..,  wu)T ,  w  =  (wi, ..,  wu)T.  Since  we  want  simplicity  and  efficiency  in 
solving  the  algebraic  equations  corresponding  to  the  implicit  part  we  will  consider  only 


(98) 

(99) 
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diagonally  implicit  Runge-Kutta  (DIRK)  schemes  for  the  source  terms  (aVJ  =  0.  for  j  >  i). 
The  use  of  a  DIRK  scheme  is  enough  to  assure  that  the  operator  Wxf  is  always  evaluated 
explicitly.  The  type  of  schemes  described  can  be  represented  with  a  compact  notation  by 
a  double  Butcher  tableau 


A 


c 

A 

T 

(J1 

where  the  coefficients  c  and  c  are  given  by  the  usual  relation 

i—  1  i 

Ci  ^  ^  CLjj ,  C{  ^  ^  CLjj .  (100) 

3= 1  J=1 

We  refer  to  Pareschi  and  Russo  (2005)  for  a  discussion  of  the  stability  requirements 
and  the  derivation  of  the  order  condition  of  IMEX  schemes  up  to  third  order.  Let  us 
remark  that  order  conditions  are  simplified  under  the  assumptions  c  =  c  and  w  =  w.  The 
first  assumption  however  prevents  the  scheme  from  being  asymptotic  preserving  unless 
the  initial  data  is  well-prepared,  namely  consistent  with  the  limiting  equilibrium  system. 

As  shown  in  Pareschi  and  Russo  (2005)  it  is  easy  to  prove  the  following: 

Lemma  4  If  all  diagonal  element  of  the  triangular  coefficient  matrix  A  that  characterize 
the  DIRK  scheme  are  non  zero,  then 

lim  F^=M[F{%  (101) 

£ — ^0 

Now  if  we  multiply  the  IMEX  method  by  the  collision  invariants  <f>{v)  =  l,v,v2  and 
integrate  the  result  in  velocity  space  we  obtain 


F®(f>(v)  dv 


fn+1cf)(v )  dv 


VxF^(f>(v)  dv 
VxF^^)(v)  dv. 


(102) 

(103) 


Thus  if  (101)  holds  true  the  higher  order  moments  of  the  F®  can  be  computed  as  function 
of  mass,  momentum  and  temperature  of  F^  and  we  get  the  explicit  Runge-Kutta  scheme 
applied  to  the  corresponding  system  of  compressible  Euler  equations.  We  can  state: 

Theorem  4  If  det  A  f  0,  in  the  limit  e  — >  0,  the  IMEX  scheme  (98)-(99)  applied  to 
system  (97)  becomes  the  explicit  RK  scheme  characterized  by  (A,w,c)  applied  to  the  limit 
Euler  system  (23). 

Note  however  that  the  above  results  do  not  guarantee  that 

lim  fn+1  =  M[fn+1]. 
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The  latter  property  is  achieved,  for  example,  if  we  assume  that 

aui  =  Wi,  aui  =  Wi,  i  =  (104) 

with  auv  0.  In  fact,  as  a  consequence  we  have  fn+1  =  F(u\  and  lim£_>0  =  M[F^}. 

As  for  the  corresponding  DIRK  methods  such  schemes  are  referred  to  as  stiffly  accurate. 

Let  us  remark  that  the  IMEX  scheme  (98)-(99)  can  be  solved  explicitly  despite  the 
nonlinearity  of  M[f],  In  fact  since  M[F^]  depends  only  on  mass,  momentum  and  tem¬ 
perature  of  the  solution,  thus  on  the  low  order  moments  of  F^\  it  can  be  evaluated 
directly  from  the  explicit  scheme  (102)-(103).  This  property  has  been  used,  for  example, 
in  Pieraccini  and  Puppo  (2007);  Filbet  and  Jin  (2010)  to  implement  efficiently  IMEX 
methods. 

Several  AP  IMEX  schemes  up  to  third  order  can  be  found  in  Pareschi  and  Russo 
(2005,  2000a).  In  Tables  9  and  10  we  report  two  examples  of  second  order  methods, 
PR(2,2,2)  method,  satisfying  w  =  w  and  det(A)  0,  and  ARS(2,2,2),  satisfying  c  =  c 
and  the  stiffly  accurate  requirement.  A  modification  of  such  schemes  in  order  to  achieve 
uniformly  accuracy  in  £  has  been  studied  in  Boscarino  and  Russo  (2009). 

4.3.2  Application  to  the  Boltzmann  equation 

In  the  sequel  we  will  describe  how  to  extend  the  IMEX  Runge-Kutta  methods  to  the  full 
Boltzmann  equation  (1).  We  follow  the  methodology  introduced  in  Filbet  and  Jin  (2010); 
Dimarco  and  Pareschi  (2011). 

Our  scope  is  to  avoid  the  implicit  solution  of  the  collision  term  of  the  Boltzmann 
equation.  To  this  goal  we  can  reformulate  the  collision  operator  adding  a  BGK  (or  another 
linear  kinetic  model  which  can  be  easily  inverted)  as  a  penalization,  exactly  as  shown  in 
Section  3.2.1.  Let  us  mention  that  a  similar  approach  has  been  previously  used  in  other 
contexts  (see,  for  example,  (Smereka,  2003))  with  the  goal  to  adopt  a  time  step  which  is 
much  larger  than  O(e)  in  regions  close  to  the  fluid  dynamic  limit. 

The  Boltzmann  equation  can  be  rewritten  as 

Btf  +  !>  ■  VJ  =  D(/)  +  -(M[f]  -  /),  (105) 

£  £ 

with 

m(f)  =  P(f,  f)  -  mMI/1,  P(f,  f)  =  Q(f,  /)  +  tif. 

Clearly  now  M[f]  is  non  constant  in  time  during  the  relaxation  process.  One  can  now  use 
a  numerical  scheme  in  which  only  the  BGK  term  M[f]  —  f  is  treated  implicitly,  while  the 
term  g(f )  describing  the  deviations  from  a  BGK  behavior  and  the  convection  term  Vx/ 
are  treated  explicitly. 

The  IMEX  scheme  now  take  the  form 

i—  1  v 

FW  =  r  +  A ( -g(F(j) )  -  V  ■  +  A tJ2  Oij-(M[FW]  -  Fij))  (106) 

3= 1  E  i=1 

V  V 

fn+ 1  =  fn  +  AtJ^Ui  (-0(FW)  -v  VXF^  +  A  tJ2^i-(M[F{i)]  -  F^).  (107) 
i= 1  '  j  I 
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Table  9:  Tableau  for  the  explicit  (left)  implicit  (right)  PR(2,2,2)  scheme 
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Table  10:  Tableau  for  the  explicit  (left)  implicit  (right)  ARS(2,2,2)  scheme 


If  we  integrate  the  above  scheme  against  <f>(v)  =  l,v,v2  we  obtain  again  system  (102)- 
(103).  It  is  a  simple  exercise  to  verify  that  Lemma  4  and  Theorem  4  apply  also  to  this 
reformulated  problem.  Thus  we  obtain  AP  method  for  the  Boltzmann  equation  that  can 
be  evaluated  explicitly.  Here,  however,  the  additional  difficulty  is  given  by  the  fact  that 
we  are  integrating  explicitly  the  stiff  term  fig{f)/e  which  may  lead  to  unstable  schemes  if 
not  properly  treated.  The  stiffly  accurate  conditions  (104)  are  essential  in  such  a  situation 
in  order  to  obtain  unconditionally  stable  AP  schemes. 

Finally  let  us  emphasize  that  the  usual  stability  requirements  may  not  suffice  to  guar¬ 
antee  nonnegativity  of  the  solution  fn+1  starting  from  a  nonnegative  initial  data  fn.  We 
refer  to  Dimarco  and  Pareschi  (2011)  for  further  details. 
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